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1.    INTRODUCTION 

The  USAF  Western  Space  and  Missile  Center  is  located  along  the 
southern  California  coast  at  Vandenberg  AFB,  where  large  Atlas  and 
Titan  rockets  are  used  to  launch  polar  orbiting  satellites.   Since 
the  rockets  are  fueled  by  tanker  trucks,  USAF  has  directed  much 
effort  toward  avoiding  or  mitigating  the  potential  for  large, 
accidental,  ground  level  releases  of  toxic  fumes. 

If  released,  such  fumes  would  disperse  downwind  according  to  the 
detailed  mechanics  of  the  release,  nature  of  the  propellants,  and 
prevailing  atmospheric  conditions  and  terrain.   Even  during  normal 
launch,  massive  amounts  of  propellants  are  released  into  the 
atmospheric  boundary  layer  via  the  rocket  exhaust  cloud.   Thus, 
USAF  focuses  part  of  its  mitigation  effort  on  field  studies  and 
computer  f orecast/nowcasts  of  such  cold  and  hot  spills  to  prepare 
toxic  hazard  corridors  (THCs)  for  base  and  civilian  personnel. 

For  example,  Battelle  Corp.  conducted  the  Mt .  Iron  dispersion  study 
from  the  three  south  base  launch  complexes,  SLC-4,  SLC-5,  and 
SLC-6,  with  252  zinc  sulfide  tracer  releases  from  Dec.  1965  to  July 
1967  (Hinds  and  Nickola,  1967).  Until  the  recent  Lompoc  Valley 
Diffusion  Study  (LVDE)  of  potential  emissions  from  the  Hypergolic 
Stockpile  and  Storage  Facility  (HSSF)  (Skupniewicz ,  et  al.  1990, 
1991),  Mt .  Iron  was  the  only  field  study  involving  tracer  releases 
conducted  from  South  Vandenberg.  Eighty-eight:  releases  were 
conducted  from  SLC-4,  25  from  SLC-5,  and  139  from  SLC-6. 

Since  SLC-6  is  presently  decommissioned,  SLC-4  is  now  the  main 
launch  pad  for  large  Vandenberg  rockets,  Thus,  in  an  effort  to 
assess  possible  replacements  for  the  currently  used  Ocean  Breeze/ 
Dry  Gulch  (OBDG)  type  regression  equations  (Haugen  and  Taylor, 
1963),  more  physically  based  computer  models  such  as  AFTOX  and 
HOTMAC/RAPTAD  (Kunkel  and  Izumi,  1991;  Yamada  and  Bunker,  1991) 
have  been  used  to  simulate  Mt .  Iron  releases  from  SLC-4.   The  Naval 
Postgraduate  School  (NPS)  has  made  available  an  eight  case  subset 
of  the  Mt.  Iron  data  set  to  the  Vandenberg  modeling  community  and 
has  made  efforts  to  coordinate  the  various  simulation  studies. 
This  report  describes  the  use  of  LINCOM/RIMPUFF,  a  dispersion  model 
from  RISO  National  Laboratory  of  Denmark,  to  simulate  the  above 
eight  representative  cases. 


2.   MOUNTAIN  IRON  DATA  DESCRIPTION 


2.1    The  Vandenberg  AFB  Domain 

Vandenberg  AFB  is  located  about  100  miles  northwest  of  Los  Angeles 
on  the  California  coast.  The  Vandenberg  domain  features  generally 
hilly  terrain,  punctuated  by  the  east-west  oriented  Santa  Ynez 
River  Valley  which  separates  North  and  South  Vandenberg.  The  South 


Vandenberg  terrain  is  generally  rougher  than  that  of  North 
Vandenberg  with  ridge  to  canyon  elevation  differences  ranging  from 
around  100  to  400m.  At  Pt .  Arguello  (near  SLC-6)  the  coastline 
slants  sharply  eastward  away  from  the  generally  north-south 
direction.  This  is  repeated  further  south  beyond  Vandenberg  by  a 
similar  coastline  turning  feature  at  Pt.  Concepcion.  Honda  Ridge 
on  South  Vandenberg  has  a  high  point  at  Tranquillon  Peak  of  over 
600m  and  an  average  height  of  about  200m.  It  represents  the 
northern  edge  of  the  coastal  Santa  Ynez  Mountains.  At  an  average 
elevation  of  about  60m,  Honda  Canyon,  roughly  2km  wide,  separates 
Honda  Ridge  from  Target  Ridge,  the  first  ridge  south  of  SLC-4  (see 
fig.  1). 

The  climatological  flow  is  from  the  northwest,  controlled  by 
subsiding  anti-cyclonic  flow  around  the  chronic  eastern  Pacific 
High.  Adiabatic  compression  of  this  subsiding  air  mass  creates  the 
typically  strong,  low  California  coastal  subsidence  inversion  seen 
at  Vandenberg.  The  northwesterly  flow  is  augmented  by  the  summer 
time  California  coast-to-Central-Valley  monsoon,  as  well  as  the 
local  Seabreeze.  However,  the  turning  coastline  at  Pts.  Arguello 
and  Concepcion  induces  a  local  divergence  of  the  climatological 
northwesterlies .  Moreover,  Honda  Ridge,  together  with  the 
typically  low  subsidence  inversion,  tends  to  divert  daytime 
northwesterly  flows  to  the  west,  so  that  winds  over  SLC-6  are 
usually  channeled  to  a  narrow  range  about  the  north-south  direction 
(see  Vandenberg  Meteorology  and  Plume  Dispersion  Handbook,  Kamada 
et  al,  1989,  hereafter  referred  to  as  the  Handbook). 

2.2    General  Data  Description 

113  tests  were  conducted  at  SLC-4  and  SLC-5  during  Phase  I  of  Mt. 
Iron  in  the  first  half  of  1966.  In  the  experimental  design, 
fluorescent  zinc  sulfide  tracer  particles  with  a  mean  solid 
diameter  of  2.5  /im  was  released  as  a  wet  aerosol  fog  from  a  test 
stand  at  a  height  of  12  ft,  usually  for  30  minutes.  Sample  tracer 
dosages  were  collected  on  membrane  filters  and  assayed  for 
fluorescence  by  alpha  emission,  using  a  Rankin  counter.  Due  to  the 
rugged  terrain,  the  sampling  lines  were  limited  to  existing  paved 
and  gravel  roads  as  shown  in  fig.  2.  Sampler  spacing  ranged  between 
approximately  160  -  640  meters  (0.1  -  0.4  miles),  depending  on 
distance  from  the  release  point.  The  160m  spacing  was  used  for 
roads  near  the  release  site.  A  320m  spacing  was  used  along  the 
eastern  section  of  Santa  Ynez,  Honda  Canyon,  Honda  Ridge,  and 
Tranquillon  Mountain  Rds.  The  640m  spacing  was  reserved  for  the 
section  of  the  Sudden  Coast  Rd.  from  Pt.  Arguello  stretching  south 
by  southeast  to  Jalama  Beach  . 

For  this  data/model  comparison,  we  selected  eight  representative 
cases.   The  cases  were  chosen  on  the  basis  of  completeness  of  the 
data  set  collected  by  Battelle  during  Mt.  Iron  and  the  significance 
of  the  presumed  flow  type  as  defined  in  the  Handbook. 
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Fig.  1    Terrain  Map  of  South  Vandenberg.   Inner  rectangle  shows 
the  100m  inner  grid  area  for  RIMPUFF.  Elevations  are  in  meters. 


VAFB   Roads    and  LIN  COM  Model  Domain 
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Fig.  2  Sampling  routes  over  South  Vandenberg  for  Mt .  Iron.  Inner 
rectangle  shows  the  100m  inner  grid  area  for  LINCOM.  Weather 
towers  are  numbered. 


The  available  Mt .  Iron  sensor  array  consisted  of  towers,  sondes, 
and  bag  samplers.   There  were:  1)  wind  speeds,  directions,  and 
direction  standard  deviations  from  19  towers  over  7.5,  30,  60,  and 
90  minute  averages  (shown  in  Tables  2.1.1,  2.1.2  and  fig.  2),   2) 
tethersonde  temperatures  at  approximately  20  meter  intervals  up  to 
260  m.  The  tethersonde  (often  referred  to  as  a  wiresonde)  was 
located  at  the  release  site,  either  VIP-1  at  the  northwest  corner 
of  SLC-4   or  Area   529  at  the  southwest  corner  of  SLC-5.   3) 
twice-a-day  NWS  rawinsondes  from  Building  22,  4)  supplemental  rawin 
or  radiosondes  from  VIP-1,  Scout  D,  upper  Honda  Canyon,  and  the 
Boathouse.   The  rawinsondes  were  released  as  often  as  every  hour 
for  3-4  hours,  usually  commencing  one  hour  before  tracer  release. 
5)   bag  sampling  lines  were  deployed  along  Mesa,   Coast,  Surf, 
Arguello,  Santa  Ynez,  Bear  Creek,  Folo,  Tank,  Target  Ridge,  Honda 
Canyon,  Tranguillon  Mt .  ,  Honda  Ridge,  and  Sudden  Coast  roads,  using 
a  total  of  "  330  bag  samplers.   The  Queen  Air  aircraft,  used  for  24 
flights  during  Phase  II,  was  deployed  only  for  case  110  of  the 
eight  cases  selected  for  study.   For  case  110  the  plume  centerlines 
derived  from  the  aircraft  data  were  essentially  identical  with 
those  from  ground  sampling. 

Table  2.1.1   Tower/Sonde  sites 


Program  Index. 

1 

2 

3 

4 

5 

6 

7 

8 

9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 

The  asterisks  denote  rawinsonde  launches  at  the  same  site  as  the 
tower 


VAFB  Site  #. 

301 

056 

014 

012 

101 

Target  Ridge 

011 

300 

054 

100 

VHF 

010 

009 

055 

057 

(Boathouse,  BH) 

* 

013 

(Honda  Canyon,  HC) 

* 

200 

(Scout-D,  SD) 

* 

VIP-1 

*  (Source) 

B22 

(Bldg.22) 

* 

Table  2.1.2   Wind  Speed,  Direction,  9,  and 
Standard  Deviation,  a0,  for  Main  Sensor  Sites 


Site 

28 

31 

48 

55 

87 

90 

91 

110 

Case 

VIP-1 

6.7 

1.5 

5.0 

5.2 

3.9 

3  .  3 

2.7 

V  m/s 

285 

330 

318 

321 

250 

325 

360 

270 

9   deg 

17 

5 

6 

8 

6 

15 

°q   deg 

VHF 

3.8 

7.4 

7.9 

4.4 

5.1 

2.8 

5.0 

V 

294 

341 

333 

275 

303 

340 

296 

9 

7 

13 

14 

11 

5 

5 

°Q 

Target 

1.2 

9.6 

2  .  3 

5.7 

1.  1 

5.8 

6.8 

2  .  0 

V 

Ridge 

269 

4 

55 

339 

271 

298 

324 

323 

9 

3 

19 

6 

11 

12 

9 

6 

9 

°e 

WT101 

4.  1 

7.4 

4  .  1 

7.4 

V 

300 

330 

325 

303 

251 

300 

315 

325 

9 

8 

15 

8 

15 

2 

12 

14 

4 

ae 

Teleme 

try 

2  .4 

4  .  1 

4.9 

2.0 

6.0 

5.7 

2.0 

V 

Peak 

321 

335 

309 

318 

310 

313 

16 

9 

5 

8 

10 

4 

12 

12 

4 

°Q 

Ion 

2.5 

2  .  0 

1.0 

2.9 

3.4 

2.9 

3.5 

4.  1 

V 

Sounder 

259 

328 

235 

288 

256 

1 

360 

249 

9 

Honda 

Cyn 

5 

4 

2 

6 

7 

6 

7 

8 

°e 

LaSall 

e 

2  .4 

1.  0 

2.9 

2.6 

4.0 

2.  6 

2.  1 

V 

Canyon 

Rd 

328 

212 

298 

268 

315 

316 

275 

9 

Honda 

Cyn 

5 

9 

5 

8 

5 

4 

ae 

WT200 

1.9 

2.  3 

0.5 

1.2 

0.8 

2.7 

V 

274 

220 

10 

256 

254 

184 

105 

259 

9 

4 

5 

7 

8 

2 

2 

5 

°e 

WT301 

2.2 

7.4 

3.9 

7.8 

7.6 

2.8 

V 

339 

339 

340 

204 

1 

15 

329 

9 

5 

15 

15 

8 

16 

15 

6 

°Q 

Boat 

3.5 

7.5 

7.0 

2.3 

13.0 

14.5 

4.6 

V 

House 

280 

350 

315 

340 

213 

309 

313 

239 

9 

7 

17 

13 

7 

26 

29 

9 

°e 

2.3    Case  Study  Selection 

As  the  Mt.  Iron  study  gained  momentum,  hourly  rawinsondes  became 
available  at  up  to  four  locations.   However,  this  was  not  true  in 
the  early  going,  and  the  wiresonde  at  the  release  site  was  not 
operative  for  the  first  two  dozen  or  so  cases.   Thus,  upper  air 
data  during  the  early  stages  of  the  study  was  rather  sparse.   Even 
some  of  the  ground  based  towers  were  not  yet  installed.   So, 
although  the  Mt .  Iron  study  began  in  the  beginning  of  December 
1965,  our  first  study  case  with  fairly  complete  data  was  case  28 
which  took  place  in  February  of  1966. 

Many  of  the  sensors  also  were  inoperative  at  various  times, 
particularly  the  wiresonde.   Since  release  site  winds  and 
temperature  structure  are  usually  the  most  critical  inputs  when 
simulating  plume  transport,  we  tried  to  limit  our  selections  to 
those  cases  where  wiresonde  data  was  available. 

Within  the  available  data,  our  other  primary  consideration  was  to 
sample  as  many  flow  conditions  as  possible.  Table  2.2.1  summarizes 
ten  significant  flow  types,  as  characterized  in  the  Handbook. 


Table  2.2.1   Surface  Flow  Types  and  Typical  Conditions 


FLOW 

FLOW 

ANGL 

TYPE 

ALOFT 
052 

AT 

1. 

Moderate 

SE  over 

190- 

South 

mixed* 

310 

westerly 

layer 

2. 

Weak 

Weak 

240- 

Westerly 

290 

HEIGHT 


Low 


Low 


USUALLY 
OCCURS 


COMMON  CAUSES 
AND  FEATURES 


spring  Seabreeze/easterly 
&  fall    diurnal  balance, 
mesoscale  eddy 


winter , 
or  with 
summer 
stratus 


Seabreeze  with 
mild  land-sea 
temperature 
difference 


3  .  Summer 
North 
westerly 


Weak 


4 .  Evening      Weak 
Northerly 

5.  Stagnation   Weak 


6.  Nocturnal 
Easterly 


290- 
340 


Low/Mid 


340 
070 

000- 
360 


-   Mid 


varies   040- 
160 


After 
noon 


Non- 
winter 


None        Near 
(surface)   dawn 


None       winter 
(Surface) 


Pacific  High  & 
Monsoon 


Veering  behind 
Seabreeze  front 

Weak  Seabreeze/ 
slope  drainage 
balance. Usually 
not  stationary 

Slope  drainage, 
winter  monsoon, 
&  diurnal 


7.  Synoptic 
Easterly 


Strong   040-   High  or     spring 
E        160    None        &  fall 


Warm,  dry  Santa 
Ana 


8.  Storm 
Southerly 

9.  Winter 
North 
westerly 


10.  Baja 

Southerly 


SW  w/  090- 

upper  200 
trough 

Strong  290- 

NW  03  0 


090- 
170 


High  or 
None 


Mid/High 
or  None 


winter 


Mid 


Frontal 
passage 


After    Postfrontal  or 
noon  or   Nevada  Low  if 
storms   strong,  Pac  High/ 
Seabreeze  if  weak 

fall     Surface  low  over 
Baja  California 


Low  =  100m 
Mid  =  400m 
High  =  600m 


March 

29 

April 

21 

June 

13 

June 

21 

June 

22 

June 

22 

With  the  above  considerations  in  mind  we  selected  the  following 
cases  from  the  1966  data  period: 


Case  Date 

28  February  16 

3  1  February  24 

48 
55 
87 
90 
91 
110 


All  the  flow  types  in  Table  2.2.1,  save  for  the  Baja  Southerly,  can 
be  viewed  within  the  context  of  seasonal  and  diurnal  variations  in 
the  Vandenberg  area  according  to  Table  2.2.2,  also  taken  from  the 
Handbook.   Perhaps  as  many  as  five  or  six  of  the  above  flow  types 
are  represented  among  the  eight  cases.   As  is  consistent  with 
Vandenberg  climatology,  the  vast  majority  of  the  Mt .  Iron  flows 
were  either  Seabreeze  or  synoptic,  post-frontal  northwester lies . 
These  are  represented  in  cases:  31,  55,  and  90.   Weaker  Seabreeze 
westerlies  are  also  present  in  cases,  28  and  110.   The  evening 
northerly  and  perhaps  partial  stagnation  is  indicated  in  cases,  48 
and  91.   The  remaining  case,  87,  is  a  Seabreeze  southwesterly, 
perhaps  initiated  by  the  positioning  of  the  edge  of  the  local 
stratus  deck  rather  than  competition  between  the  Seabreeze  and  a 
synoptic  easterly. 

Storm  southerlies,  nocturnal  easterly  drainage,  synoptic 
easterlies,  and  Baja  Southerlies  are  not  indicated  in  the  available 
data.   The  first  of  the  three  cases  which  might  be  construed  as 
storm  southerlies  does  not  appear  in  the  Mt.  Iron  data  set  until 
MI-175  from  SLC-6  in  January  1967  during  Phase  II.  Since  data 
collection  during  storm  periods  is  difficult  or  impossible,  perhaps 
this  was   intentionally  avoided  until   a  high   level   of   field 
competence  was  achieved.   The  scarcity  of  easterly  drainage  flow 
cases  during  both  Phase  I  and  Phase  II  is  not  surprising,  since 
most  data  were  taken  during  daylight  hours. 

The  upper  air  data  was  also  valuable  for  visualizing  the  temporal 
and  spatial  evolution  of  the  flow  over  Vandenberg  (e.g,  sea-breeze 
front  location,  marine  layer  depth,  drainage,  fog  or  stratus 
effects)  and  aided  in  determining  the  flow  type  analysis. 

Since  a  number  of  rarer  flow  types  characterized  in  the  Handbook 
were  not  available  within  the  data  set,  we  also  tried  to  skew  the 
sampling  to  span  all  three  available  seasons  and  some  night  time 
cases.   Thus,  cases  31,  48,  90,  and  91  were  nocturnal  releases, 


while  cases  28,  55,  87,  and  110  were  in  daytime.  Cases  28  and  31 
occurred  during  late  winter,  48  and  55  in  spring,  and  87,  90,  91, 
and  110  occurred  in  early  summer.  Tracer  was  released  from  SLC-6 
exclusively  during  the  fall  season. 

Seven  of  the  eight  cases  were  releases  from  SLC-4 ,  with  only  case 
55  being  a  release  from  SLC-5.  Though  not  from  SLC-4,  case  55  was 
chosen  because  it  was  the  only  case  indicating  strong  convection. 
In  Phase  I  there  were  some  instances  of  plumes  from  SLC-5  largely 
confined  to  flow  up  Honda  canyon  during  the  usual  mid-morning 
Seabreeze  westerlies.  But  this  pattern  was  not  seen  from  SLC-4  and 
we  did  not  include  them. 


2.4  Case  Flow  Analysis 

The  eight  Mountain  Iron  cases  extend  from  February  to  June,  with 
marine  boundary  layer  heights  varying  between  160  and  1300  meters. 
The  following  are  capsule  descriptions  of  the  meteorology. 


MI-28 

Date 
Release  time 

site 

duration 
Flow  Type 
Inversion 


16  February  1966 

1200   local 

VIP  #1 

15  min 

Seabreeze  weak  westerly  (SWW) 

Deep,  weak,  at  1300  meters? 


Flow  analysis  :  The  small  horizontal  variation  in  the 
temperature  profile  suggests  modestly  higher  off-coast  pressures 
indicative  of  a  typical  noontime,  weak  winter  Seabreeze.  Wind 
vectors  at  WT301  and  the  Boathouse  suggest  that  the  flow  diverges 
around  the  high  terrain  across  Honda  Ridge,  where  hill  induced 
pressure  gradients  create  the  surface  layer  acceleration,  evident 
at  the  500m  level  at  WT055.  The  veering  indicated  at  Telemetry 
Peak  and  WT014  are  probably  due  to  local  divergence  initiated  by 
the  high  terrain.  Even  at  these  modest  wind  speeds  the  evident 
flow  distortion  generated  by  Honda  Ridge  suggests  a  lower,  stronger 
inversion  than  is  indicated  by  sounding  records. 


MI-31 


Date,  release  time 

Release  time 
site 
duration 

Flow  Type 

Inversion 


24  February  1966 

1845  local 

VIP  #1 

15  min 

Synoptic  Northwesterly  (NW) 

Unknown,  probably  deep 


10 


Flow  analysis  :  Wiresonde  temperatures  were  available  only  to 
220m.  Hence,  the  lack  of  upper  wind  or  cloud  data  precludes 
certainty.  However,  the  strong  northwesterly  flow,  wind  speeds  and 
time  of  year  suggest  post-frontal  forcing  with  modest  veering 
within  Honda  Canyon  and  surface  layer  acceleration  across  Honda 
Ridge.  The  lack  of  flow  divergence  to  the  westward,  lower  portion 
of  Honda  Ridge  also  suggests  a  high  inversion,  consistent  with 
post-frontal  conditions.  The  abscence  of  backing  once  the  flow  has 
passed  beyond  the  ridge  is  consistent  with  this  interpretation. 


MI-48 


Date 
Release  time 

site 

duration 
Flow  Type 


29  March  1966 

2315  local 

VIP  #1 

5  min 

Nocturnal  Drainage  (ND) 


Inversion  :   Elevated  subsidence  type  at  400m  near  release 

site  but  with  large  variability.   For  example  the  inversion  above 
the  boathouse  appears  low,  perhaps  due  to  radiation  fog. 

Flow  analysis  :  Tower  winds  show  veered,  northwesterly 
remnants  of  a  Seabreeze  disturbed  by  more  veering  across  Honda 
Ridge.  The  inversion  lid  forces  the  flow  '/est  over  the  lower 
portion  of  Honda  Ridge.  South  of  the  ridge,  continuity  forces  the 
winds  east  again,  as  seen  at  the  Boathouse.  Downslope  drainage 
flow  into  the  canyons  is  seen  at  low  levels,  with  complicated 
reverse  flow  due  to  radiation  fog  which  forms  in  the  larger 
drainage  pool  areas  such  as  upper  Honda  Canyon. 


MI-55 

Date 
Release  time 

site 

duration 
Flow  Type 

Inversion 


21  April  1966 
1109  local 
Bldg  529  (SCOUT) 
30  min 
Northwesterly 

Weak  subsidence  type  at  650m 


Flow  analysis  :  Well-mixed  boundary  layer  up  to  the  inversion 
with  winds  maximizing  at  200-300m,  then  decreasing  aloft.  Humidity 
profiles  indicate  cloud  cover  over  Honda  Canyon  and  Scout  D.  In 
view  of  the  unstable  lapse  rates  indicated  by  both  the  rawin  and 
wiresondes,  this  suggests  spring  cumulus  convection,  consistent 
with  the  short  plume  footprint.  Flow  distortion  across  Honda 
Ridge,  similar  to  MI-48,  is  still  evident  but  weaker,  perhaps  due 
to  the  higher,  weaker  inversion. 
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MI-87 

Date 
Release  time 

site 

duration 
Flow  Type 

Inversion 


13  June  1966 

1310  local 

VIP  #1 

30  min 

Southwesterly  Seabreeze 

Low  subsidence  type  at  250m. 


Flow  analysis  :  The  flow  is  weakened  and  backed  from  its  usual 
summertime  pattern  near  the  surface.  However,  rawinsondes  and 
towers  along  the  ridge  lines  show  that  the  higher  level  winds 
retained  the  usual  northwesterly  Seabreeze  character  for  this  time 
of  day.  This  is  the  opposite  pattern  from  that  expected  for 
southwesterlies  induced  by  inland  high  pressures,  the  kind  which 
give  rise  to  hot  e  terly  Santa  Anas  further  south.  In  this  case 
the  Seabreeze  has  ought  a  fog  front  with  it,  leading  to  a  low 
level,  thermally  based,  southwesterly  "fog  breeze". 


MI-90 

Date 
Release  time 

site 

duration 
Flow  Type 

Inversion 


21  June  1966 

2300  local 

VIP  #1 

30  min 

Northwesterly  summer  Seabreeze 

High  subsidence  type  at  600m. 


Flow  analysis  :  The  high  inversion,  clear  skies  and  strong 
winds  this  late  into  the  evening  suggest  some  frontal  augmentation 
of  the  usual  Seabreeze.  The  usual  flow  distortion  over  and  around 
Honda  Ridge  is  again  shown  by  the  wind  vectors  at  WT301,  055,  and 
057.  Tripling  of  the  speed  at  WT055  suggests  strong  flow  funneling 
near  the  inversion,  as  well  as  acceleration  over  hills,  due  to  hill 
induced  pressure  gradients  (Jackson  and  Hunt,  1979)  .  The  anomalous 
southerly  flow  at  the  mouth  of  Honda  Canyon  is  hard  to  explain, 
except  as  instrument  error.  However,  we  speculate  that  the  high 
winds  at  the  Boathouse  may  be  due  to  a  lowered  inversion  height  and 
conseguent  high  speed  funneling  south  of  Honda  Ridge.  The  boundary 
layer  may  be  lower  south  of  the  ridge  because  the  ridge  obstructs 
and  dams  the  northwesterly  flow,  inducing  a  Scorer  type  of 
hydraulic  jump  across  the  ridge.  Rippling  of  the  inversion  base  by 
leeside  gravity  waves  may  augment  this  effect. 
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MI-91 

Date 
Release  time 

site 

duration 
Flow  Type 

Inversion 


22  June  1966 

0203  local 

VIP  #1 

30  min 

Evening  Northerly 

Mid-level  subsidence  type  at  400m 


Flow  analysis  :  This  release  occurred  three  hours  after  case 
90.  In  the  interim  winds  have  veered  from  northwesterly  to 
northerly  and  weakened  somewhat,  due  to  the  lack  of  thermal 
forcing.  This  follows  the  summer  Seabreeze  pattern  outlined  in  the 
Handbook.  According  to  rawinsondes,  the  veering  occurred  around 
0000  local  time  and  appeared  more  evident  above  600m,  while  the 
inversion  base  fell  sharply  (by  ~250m) .  These  changes  may  be 
consistent  with  post-frontal  relaxation  of  the  hydraulic  jump 
mentioned  in  the  case  90  analysis.  Even  so,  northerly  forcing 
remained  intense  enough  to  avoid  any  drainage  or  higher  level 
easterly  return  flow  tendencies,  as  often  appear  later  at  night. 
Again,  the  usual  flow  distortion  over  and  around  Honda  Ridge  shows 
up  in  the  wind  vectors  at  WT301,  055,  and  057.  And  the  wind  vane 
at  the  mouth  of  Honda  Canyon  still  seems  in  error. 


MI-110 

Date 
Release  time 

site 

duration 
Flow  Type 

Inversion 


22  June  1966 

1055  local 

VIP  #1 

30  min 

Seabreeze  Weak  Westerly 

strong  subsidence  type  at  160m 


Flow  analysis      :   This  release  occurred  eight  hours  after  case 
91.    In  the  interim  winds  have  veered  completely  around  to 
westerly.   Inland  heating  has  burnt  back  the  overnight  stratus 
cover  to  near  shore,  above  which  remains  the  ever  present 
subsidence  inversion.   However,  the  Seabreeze  is  still  weak,  due  to 
the  limited  heating  period,  nor  has  Coriolis  induced  veering 
really  begun  yet.    This  follows  the  summer  Seabreeze  pattern 
outlined  in  the  Handbook.   According  to  rawinsondes,  winds  backed 
rapidly  above  the  inversion  to  a  light  southeasterly,  also  typical 
of  the  general  Seabreeze  circulation  pattern  at  Vandenberg. 
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3.    LINCOM-RIMPUFF  MODEL  DESCRIPTION 


3 . 1   LINCOM 


There  are  two  general  classes  of  models,  diagnostic  (steady  state) , 
and  prognostic  (time  marching) .  Diagnostic  codes  are  suited  for 
nowcasts  where  immediate  toxic  hazard  corridors  (THCs)  are  needed 
for  accidental  spills.  Prognostic  codes  are  more  suited  for 
forecasts  of  hypothetical  spills  or  planned  launches.  This 
distinction  blurs,  if  a  diagnostic  model  is  fast  enough  to  provide 
freguent  updates  based  on  perhaps  a  five  minute  or  better  basis,  or 
if  f aster-than-real-time  prognostic  codes  are  run  continuously  and 
frequently  nudged  by  new  input  data.  LINCOM  is  a  five  level, 
diagnostic  boundary  layer  flow  model,  which  when  combined  with  a 
puff  dispersion  model  such  as  RIMPUFF,  is  fast  enough  to  allow  such 
updates  on  a  two  minute  basis.  LINCOM  1.0  itself  requires  about 
one  minute  for  a  typical  case  at  500m  resolution  over  a  40  x  60  km 
domain  on  a  386/33  PC. 

LINCOM  was  developed  by  lb  Troen  at  the  RIS0  National  Laboratory  in 
Roskilde,  Denmark  (Troen,  1986).  Other  recent  contributers  have 
been  George  Lai,  Ray  Kamada,  and  Torben  Mikkelsen. 

LINCOM  uses  linearized  versions  of  the  u,  v,  and  w  components  of 
the  momentum  equation,  the  boundary  layer  mixing  form  of  the  mass 
continuity  equation  and  the  equation  of  state.  LINCOM  versions  1.0 
and  1.1  both  neglect  the  temperature  equation  (also  termed  the 
thermodynamic  energy  equation).  Thus,  neither  stable  nor  unstable 
conditions  are  treated  within  the  governing  equations.  However, 
through  an  objective  analysis  scheme  applied  after  the  solutions 
from  the  governing  equations  are  obtained,  LINCOM  is  designed  to 
match  the  tower  data  exactly,  with  dynamic  interpolation  between 
towers.  That  is,  at  grid  points  between  the  tower  inputs,  LINCOM 
adjusts  the  winds  by  combining  linearized  dynamics  with  a  terrain 
modified,  inverse-distance-square  based  objective  analysis.  So 
LINCOM  conforms  to  the  actual  winds  and  the  thermodynamic  and 
dynamic  forcings,  to  the  level  of  accuracy  and  resolution  provided 
by  the  towers.  Since  LINCOM  operates  at  five  levels  within  the 
boundary  layer  (default  values  are:  6',  54',  0.2  Zi,  0.5  Zi,  and 
0.8  Zi)  ,  levels  above  the  towers  use  similarity  based 
extrapolations. 

LINCOM  1.0  was  used  for  the  Mt .  Iron  study.   However,  in  version 
1.1,  the  resultant  wind  field  is  not  completely  anchored  to  the 
tower  data,  as  is  discussed  below. 

LINCOM  is  speed  optimized  by  two  techniques.   First,  we  transform 
the  Vandenberg  terrain  from  grid  point  to  Fourier  spectral  space 
wherein  the  governing  eguations  are  solved.   In  this  way  the  same 
resolution  can  be  achieved  using  far  fewer  sine/cosine  waves  as 
grid  points.   Second,  instead  of  numerically  solving  the  governing 
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equations  at  each  grid  point,  linearization  lets  us  pre-calculate 
symbolic  solutions  so  that  only  the  coefficients  for  each  wave 
number  need  further  adjustment.  Together  these  two  techniques  allow 
more  than  an  order  of  magnitude  speed  increase  over  standard  finite 
difference  methods. 

In  addition,  for  linearized  neutral  flow,  Lai  (unpublished)  showed 
recently  that  the  total  flow  can  always  be  expressed  as  a  linear 
combination  of  independent  orthogonal  components  of  the  mean  and 
perturbed  flows  along  just  the  x  and  y  axes.   Since  the  orthogonal 
flows  are  pre-calculable,  we  can  simply  find  the  mean  wind  speed 
and  direction  which  best  matches  the  Vandenberg  tower  winds  and 
apply  them  to  pre-calculated  solutions.    In  version  1.0,   a 
"look-up"  mode  employs  a  set  of  72  wind  fields  spaced  in  three 
degree  increments  around  360°  of  arc  which  are  interpolated  to 
obtain  the  final  field.   However,  in  version  1.1,  the  best  linear 
combination  of  two  pre-calculated  orthogonal  fields  is  provided. 

Trade-offs  are  involved  in  formulating  such  a  model.   Since  LINCOM 
does  not  account  for  time  tendencies  in  the  velocity, 
dynamical  forcing  due  to  terrain  is  effectively  propagated 
instantaneously  over  the  entire  domain,  rather  than  at  some 
realistic  set  of  phase  speeds.   Linearization  of  the  governing 
equations  also  means  that  second  and  higher  order  perturbation 
terms  are  neglected.   In  a  prognostic  scheme,  such  solutions 
would  gradually  diverge  from  reality.   However,  for  a  diagnostic 
model,  where  in  effect  only  the  present  time  step  is  being 
considered,  this  issue  is  not  very  significant. 

Another  aspect  of  linearization  in  LINCOM  is  that  for  each  grid 
point,  the  vertical  velocity  is  basically  assumed  equal  to  the 
horizontal  velocity  at  that  site  times  the  sine  of  the  terrain 
slope.   This  approximation  only  holds  well  for  slopes  of  less  than 
twenty  degrees.   Thus,  the  aspect  ratio  of  the  terrain  (typical 
height/typical  length)  must  be  relatively  small  compared  to  unity. 
ACTA ' s  analysis  of  the  Vandenberg  terrain  concludes  that  at  500 
meter  resolution,  hardly  any  of  the  terrain  shows  slopes  greater 
than  twenty  degrees  (Conley  et  al.,  1990).   Thus,  the  quantitative 
limit  for  the  vertical  velocity  assumption  is  seldom  exceeded  at 
Vandenberg. 

Transformation  of  the  domain  to  spectral  space  involves  the 
assumption  of  periodic  boundary  conditions.   So  wave  energy 
leaving  one  boundary  will  immediately  re-enter  the  domain  from  the 
opposite  boundary.   For  this  reason,  a  10  km  buffer  zone  was 
created  around  the  Vandenberg  domain  in  which  the  terrain  is 
gradually  relaxed  on  all  sides  to  sea  level.   This  has  the  effect 
of  greatly  damping  such  wave  motions.    However,  the  effect  of 
periodic  boundary  conditions  cannot  be  avoided  completely  and  we  do 
not  have  a  quantitative  assessment  of  its  influence  on  the  flow  as 
yet. 
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Again,  the  effect  of  the  objective  analysis  in  LINCOM  1.0  is  to 
anchor  the  resultant  wind  field  to  the  wind  vectors  given  by  the 
Vandenberg  base  tower  system.  In  LINCOM  1.0  Monin-Obukhov 
similarity  theory  is  used  to  extend  the  tower  vectors  upward  into 
the  outer  boundary  layer.  In  effect  this  turns  LINCOM  into  a 
dynamically  based  interpolation/extrapolation  scheme  based  on  the 
tower  winds.  However,  the  tower  winds  may  be  influenced  by  effects 
which  are  more  local  than  can  be  resolved  by  the  implemented  500 
meter  grid  spacing,  such  as  buildings  and  local  terrain.  This 
constitutes  a  type  of  aliasing.  Since  the  objective  analysis  is 
not  inherently  mass  conserving,  such  local  streamline  divergence 
indicated  by  the  towers  can  cause  departures  from  the  continuity 
constraint  implemented  in  the  model  physics.  However,  since  LINCOM 
1.1  is  a  best  fit  to,  rather  than  an  anchoring  by  the  tower 
vectors,  strict  mass  conservation  is  retained,  while  avoiding  the 
aliasing  issue. 

Most  of  the  available  towers  are  clustered  near  the  coast  on 
Vandenberg  property.  Other  towers  are  available  from  the  Santa 
Barbara  Air  Quality  Management  District's  set  of  towers,  some  of 
which  are  actually  maintained  by  the  regional  oil  refineries. 
There  are  also  two  local  buoys  maintained  by  NOAA  and  towers  atop 
a  few  of  the  local  oil  drilling  platforms.  However,  these  sources 
have  not  been  incorporated  as  yet  into  the  on-line  data  stream. 


3.2    RIMPUFF 

There  are  four  classes  of  plume  diffusion  models  which  are  being 
considered  for  use  at  Vandenberg.   In  ascending  order  of 
computational  demands  they  are:  1)  steady  state  plume  model  with 
gaussian  lateral  concentration  profile,  2)  mixed  layer  scaling 
model,  3)  lagrangian  puff  model  with  gaussian  radial  concentration 
profile,   4)   lagrangian  particle  model  with  particle  motions 
governed  by  a  mean  wind  plus  a  random  turbulent  component  with  a 
gaussian  velocity  profile.   RIMPUFF  is  a  lagrangian  puff  model. 

RIMPUFF  is  structured  to  handle  multiple  simultaneous  sources. 
Release  points  can  be  located  anywhere  within  the  3-D  grid  and  can 
be  specified  individual  release  rates,  release  times  and  heat 
production.   A  series  of  puffs  is  released  to  simulate  a 
continuous  plume.   At  each  time  step,  a  book-keeping  algorithm 
tracks  the  advection,  diffusion,  and  fractional  deposition  of 
each  puff  in  accordance  with  local  meteorological  parameters. 
The  model  calculates  the  concentration  at  each  grid  point  by 
summing  the  contributions  from  all  the  surrounding  puffs.  These 
concentrations  can  simply  be  updated  or  also  accumulated  to  provide 
dosage.  The  model  output  consists  of  individual  puff  locations  and 
grid  concentrations  at  user  specified  time  intervals.   From  these 
results,  a  graphics  program  produces  puff  plots  and  concentration 
or  dose  isopleths. 
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RIMPUFF  employs  the  LINCOM  mean  flow  field  to  advect  contaminant 
puffs  downstream.  Within-puff  growth  is  controlled  by  local 
turbulence  level,  using  two-particle  relative  diffusion  theory. 
Local  turbulence  intensity  depends  on  tower  based  climatological/ 
directional  parameterizations  of  ae,  the  standard  deviation  of  wind 
direction,  for  lateral  turbulence,  and  Pasquill-Gif ford  stability 
classes  modified  for  complex  terrain  for  the  vertical  turbulence. 

RIMPUFF  treats  plume  bifurcation  in  complex  terrain  by  letting 
puffs  split  when  they  exceed  the  size  of  the  grid  spacing.   I.e., 
for  Vandenberg  an  initial  100m  puff  will  grow  to  the  500m  LINCOM 
grid  spacing  before  it  splits  horizontally  into  five  250m  puffs 
(pentafurcation) .   Vertical  trifurcation  can  occur  independently. 
The  practical  splitting  limit  on  PC  based  computers  is  several 
hundred  puff  progeny.  Mass,  momentum,  and  center  concentration  (of 
the  mother  puff)  are  conserved  in  the  splitting  process.   If  the 
mean  flow  is  not  frequently  updated  and  remains  static,  RIMPUFF 
adds  stochastic  (Langevin  based)  advection  to  simulate  puff  meander 
due  to  sub-grid  turbulence.   This  is  particularly  important  when 
puff  siblings  are  closely  spaced  shortly  after  splitting. 

Since  the  series  of  puffs  is  meant  to  simulate  a  continuous  plume, 
the  effective  stack  height  after  plume  rise  is  taken  from  standard 
plume  rise  models  for  hot  spill  cases.  100%  reflection  is  assumed 
at  a  user  specified  inversion  height  which  parallels  the  terrain. 
The  final  rise  height  for  each  puff  is  a  function  of  the 
atmospheric  stability  and  windspeed  at  the  time  and  height  of  the 
release  as  given  by  LINCOM. 

For  Mt.  Iron  the  ground  level  reflection  coefficient  was  set  to 
100%.   Dry  deposition  is  calculated  using  the  source  depletion 
concept  according  to  stability  and  windspeed.   Wet  deposition 
accounts  for  rain  intensity  and  duration  in  time  and  space,  using 
a  1/r2  objective  analysis,  where  r  is  the  distance  from  the  grid 
point  to  the  measurement  station. 

As  a  puff  model  RIMPUFF  has  advantages  over  gaussian  plume  and 
similarity  scaling  models.   It  can  handle  non-homogeneous  terrain 
during  non-stationary  conditions.   As  plumes  extend  to  farther 
distances  and  longer  travel  times,  they  are  more  likely  to 
encounter  and  respond  to  such  variations.   In  complex  terrain, 
non-gaussian  profiles  can  occur  which  RIMPUFF  can  treat,  such  as 
bifurcated   plumes,   pooling,   non-homogeneous   channeling,   and 
fanning. 

Gaussian  plume  and  mixed  layer  scaling  models  must  also  assume  some 
kind  of  profile  for  the  vertical  concentration  distribution.   A 
gaussian  vertical  profile  may  be  modified  to  account  for  surface 
reflection.   This  is  appropriate  for  short  distances  before  a 
significant  amount  of  reflection  occurs  from  an  upper  level 
inversion.   At  long  distances  a  uniform  vertical  concentration  may 
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be  assumed.   In  fact,  neither  of  these  simple  assumptions  is 
entirely  reasonable,  since  plumes  will  spread  vertically  until  an 
inversion  is  reached,  then  reflect  and  mix  downward  in  a  gradual 
and  complex  fashion.  Non-stationary  lagrangian  puffs  and  particles 
should  be  able  to  model  such  complex  phenomena  more  accurately. 

Puff  model  accuracy  is  mitigated,  especially  near  the  source,  by 
finite  puff  size  and  assumed  radially  gaussian  concentration 
profiles.  Lagrangian  particle  models  improve  on  these  aspects.  But 
they  require  a  much  larger  particle  swarm  and  more  computer  time. 
RIMPUFF's  stochastic  puff  advection  scheme  is  also  not  well  tested 
yet.    On  the  other  hand,  lagrangian  particle  models  assume  a 
gaussian  velocity  component  to  diffuse  the  particles.   Yet,  the 
vertical  velocity  distribution  has  positive  skewness  and  therefore 
is  not  gaussian  under  convective  conditions.    Baerentsen  and 
Berkowicz  (1984)  and  Misra  (1982) ,  however,  claim  better  results 
with   an   adjustable   double   gaussian   distribution   for   their 
lagrangian  models. 

4.0  DATA/MODEL  COMPARISON  METHODS 

4.1  Initialization 

With  a  model  to  model  comparison,  it  is  trivial  to  create  output 
grids  with  compatible  domain  size  and  grid  spacing  to 
quantitatively  compare  the  results.  However,  a  model  to  data 
comparison  such  as  Mt.  Iron  requires  some  manipulation  before 
quantitative  results  may  be  obtained. 

Among  the  products  of  the  Mt.  Iron  Experiment  was  a  set  of  hand- 
drawn  isopleths  defining  measured  dosages  for  each  case,  given  in 
terms  of  time-integrated  concentrations  normalized  by  release  rate 
(units  of  109sec/m  ).  However,  these  results  were  not  immediately 
compatible  with  the  format  required  for  comparison  with  LINCOM/ 
RIMPUFF  which  defines  dosages  at  distinct  points  on  a  grid.  Thus, 
the  eight  cases  from  Mt.  Iron  were  digitized  onto  a  119  x  84  evenly 
spaced  grid  at  100m  mesh  spacing,  and  formatted  for  the  commercial 
SURFER  PC  plotting  package  (Golden  Software,  1989).  This  domain 
size  was  large  enough  to  encompass  the  Mt.  Iron  and  LINCOM/RIMPUFF 
plumes  yet  small  enough  to  comply  with  the  grid  size  constraints  of 
SURFER  (see  figs.  1  and  2) . 

Inspection  revealed  that  some  of  the  hand-drawn  dose  isopleths  from 
Mt.  Iron  were  extended  beyond  the  scope  of  the  bag  sampling 
network,  presumably  to  ensure  closure  of  the  contour  lines.  Such 
extrapolation  beyond  the  measurements  in  cases:  28,  48,  87,  90,  and 
110  cannot  really  be  justified,  since  it  relies  on  subjective 
judgement.  In  these  cases  both  the  modeled  and  measured  isopleths 
were  truncated  along  the  last  measurement  line.   For  instance,  in 
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case  28,  the  measured  plume  passed  over  the  sampling  line  along 
Arguello  Road,  but  did  not  reach  Santa  Ynez  Road.  Since  the  exact 
loci  for  dose  isopleths  between  these  two  sampling  lines  was  not 
known,  the  isopleths  were  truncated  at  Arguello  Road  for  the  DICA 
and  other  comparisons. 

For  the  following  cases  the  cut-off  lines  were: 

Mt.  Iron  Case        Plume  Cutoff 


28 
48 
87 
90 
110 


Arguello 

Honda  Ridge 

Arguello  /  Santa  Ynez 

Tranquillon  /  Honda  Ridge 

Arguello 


Measured  dosages  along  the  Coast  Road  in  the  vicinity  of  Sudden 
Ranch  for  case  91  were  sufficient  to  resolve  the  plume,  whereas  for 
case  48,  measured  dosages  were  too  small  and  the  plume  was 
truncated  at  Honda  Ridge. 

4.2   DICA  and  Other  Merit  Scoring  Methods 

The  Dose  Isopleth  Correlation  Area  (DICA)  scoring  method  is  a 
conceptually  simple  performance  measure  we  created  to  compare  the 
coincidence  between  LINCOM/RIMPUFF  dose  isopleths  and  the  hand 
drawn  ones,  based  on  Mt .  Iron.  From  simple  set  theory,  for  each 
case  and  each  dose  isopleth  within  a  case,  DICA  divides  the  total 
area  of  intersection  between  each  modeled  and  measured  isopleth  by 
the  total  area  of  union.  However,  before  being  summed,  it  weights 
each  intersect  by  the  measured  to  modeled  isopleth  value  ratio. 
Thus,  if  there  are  k  number  of  isopleths  with  modeled  and  measured 
doses,  aj,  a:,  and  modeled  and  measured  areas,  Aj,  A]t  the  DICA 
algorithm  is 


k    k  k  k 

E     E  (a  /a  )A  n  A   +   E  E   (a  /a  )A  n  A 

j=l   i=l  i   j   i    j    i=l  j=l    j   i   j    i 

a  <a  i=j 


i   D 


a  <a 

j   i 


k        k 
E   A   U   E   A 
i=l   i    j=l   j 


(4.2.1) 
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This  yields  a  single  value,  ranging  between  zero  and  one,  for  the 
amount  of  overlap  between  modeled  and  measured  results.   Perfect 
congruence  yields  unity,  while  no  overlap  yields  zero. 

DICA  presents  a  fairly  stiff  test  of  "goodness  of  fit".  That  is, 
the  congruence  must  be  guite  good  to  obtain  DICA  values  above  0.5. 
Thus,  we  suggest  the  following  interpretation  of  DICA  scores, 


DICA  VALUE       FITTING  ACCURACY 


0.5   -  1 

excellent 

0.2   -0.5 

good 

0. 05  -  0.2 

fair 

0.00  -  0.05 

poor . 

In  addition  to  DICA  two  other  performance  measures  were  utilized  as 
suggested  by  Hanna  (1991),   fractional  bias, 


FB   =  (D0  -  Dm)/(D0  +  D,J    ,  (4.2.2) 

and  normalized  mean  sguare  error, 


NMSE  =    (D0  -  Dm)/D0Dm   ,  (4.2.3) 

where  D0  and  Dm  are  the  modeled  and  measured  dosages  at  each  grid 
point.   Programs  were  written  to  determine  DICA,  FB,  and  NMSE 
values  for  each  of  the  eight  cases. 

FB  is  useful  for  showing  which  plume  predicts  the  largest  dosages. 
A  zero  FB  indicates  that  the  model  predicts  on  average  the  same 
dosage  as  was  measured.  The  NMSE  should  be  behaviorally  similar  to 
DICA,  except  that  lower  NMSE  should  imply  better  model  predictions, 
while  higher  DICA  values  indicate  a  better  fit  (see  discussion  of 
relative  merits  in  section  5.2). 

All  values  which  exceeded  the  lowest  isopleth  value  of  20  are 
included.  If  either  the  RIMPUFF  or  Mt.  Iron  dosage  value  at  a 
given  grid  point  exceeded  the  minimum  threshold  value,  then  that 
grid  point  was  included  in  the  averaging  process  for  the  FB  and  the 
NMSE  scores. 

We  also  show  a  standard  comparison  of  observed  versus  modeled 
centerline  dose  exposure  levels,  including  factor  of  two  and  factor 
of  four  error  limits  (fig.  27). 
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5.   DATA/MODEL  COMPARISON  RESULTS 

5.1   LINCOM  Winds,  RIMPUFF/Data  Isopleths,  DICA  and  Other  Results 

Results  of  each  Mt .  Iron  case  and  LINCOM/RIMPUFF  simulation  are 
plotted  in  figs.  3-27.  Dose  isopleths  are  given  in  terms  of  time- 
integrated  concentrations  normalized  by  release  rate  (109  sec/m3)  . 
Locations  are  given  in  Universal  Transverse  Meridian  System  (UTMS) 
coordinates  at  100  m  resolution.  The  release  location  for  case  55 
was  at  the  mouth  of  Honda  Canyon  near  SLC-5.  The  other  seven 
releases  were  from  VIP  #1  at  the  northwest  corner  of  the  SLC-4  area 
along  the  Coast  Road. 

The  DICA,  FB,  and  NMSE  values  are  presented  in  Table  5.1  for  the 
eight  cases. 


Table  5.1 

.  1   Com 

parison 

of  model  p 

erformance  measures. 

Case 

Number 

of 

points 

D0 

dp 

(Do-Dp)2 

FB 

NMSE 

DICA 

28 

575 

420 

684 

1, 175,025 

-0.477 

4  .  082 

0.  360 

31 

1038 

231 

188 

125,231 

+0.207 

2  .  863 

0.266 

48 

1499 

178 

225 

135,991 

-0.233 

3.391 

0.268 

55 

69  2 

132 

142 

133, 124 

-0.069 

7.063 

0.  121 

87 

692 

264 

355 

265, 698 

-0.294 

2.825 

0.535 

90 

1589 

285 

460 

837,967 

-0.470 

6.375 

0.  328 

91 

2209 

177 

187 

387,625 

-0.054 

11.700 

0.369 

110 

517 

303 

470 

484,409 

-0.432 

3.395 

0.472 
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Fig.  3    Mt.  Iron  Case  28  plume  isopleths  and  tower  wind  vectors 
from  the  Mt .  Iron  Report  Phase  I  Vol.  I. 
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39.5- 


Scale 

—  4.70  m/s 


Fiq   4     LIHCOM  Case  28  (weak  Seabreeze  westerly)  wind  vectors 
shown  at  every  fourth  grid  point  (2km).   Domain  is  25  x  40  km. 
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LINCOM  vs.    ML    Iron   Case   28 
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Fig.  5    LINCOM/RIMPUFF  and  Mt.  Iron  dose  isopleths  for  Case  28  at 
100m  resolution.   Domain  is  8.4  x  11.9  km. 
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Fig.  6    Mt.  Iron  Case  31  plume  isopleths  and  tower  wind  vectors 
from  the  Mt.  Iron  Report  Phase  I  Vol.  I. 
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Fig.  7   LINCOM  Case  31  (synoptic  northwesterly)  wind  vectors  shown 
at  every  fourth  grid  point  (2km) .   Domain  is  25  x  40  km. 
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LIN  COM  vs.    ML    Iron   Case    31 
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Fig.  8    LINCOM/RIMPUFF  and  Mt.  Iron  dose  isopleths  for  Case  31  at 
100m  resolution.   Domain  is  8.4  x  11.9  km. 
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Fig.  9    Mt.  Iron  Case  48  plume  isopleths  and  tower  wind  vectors 
from  the  Mt.  Iron  Report  Phase  I  Vol.  I. 
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Fig.  10   LINCOM  Case  48  (nocturnal  drainage)  wind  vectors  shown  at 
every  fourth  grid  point  (2km).   Domain  is  25  x  40  km. 
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LINCOM  vs.    Mt.    Iron   Case    48 
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Fig.  11   LINCOM/RIMPUFF  and  Mt .  Iron  dose  isopleths  for  Case  48  at 
100m  resolution.   Domain  is  8.4  x  11.9  km. 
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Fig.  12   Mt.  Iron  Case  55  plume  isopleths  and  tower  wind  vectors 
from  the  Mt .  Iron  Report  Phase  I  Vol.  I. 
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Fig.  13  LINCOM  Case  55  (convective  Seabreeze  northwesterly)  wind 
vectors  shown  at  every  fourth  grid  point  (2km) .  Domain  is  25  x  40 
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Fig.  14   LINCOM/RIMPUFF  and  Mt .  Iron  dose  isopleths  for  Case  55  at 
100m  resolution.   Domain  is  8.4  x  11.9  km. 
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Fig.  15   Mt.  Iron  Case  87  plume  isopleths  and  tower  wind  vectors 
from  the  Mt .  Iron  Report  Phase  I  Vol.  I. 
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35 


LIN  COM  vs.    ML    Iron   Case    87 
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Fig.  17   LINCOM/RIMPUFF  and  Mt.  Iron  dose  isopleths  for  Case  87  at 
100m  resolution.   Domain  is  8.4  x  11.9  km. 
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Fig-  18   Mt.  Iron  Case  90  plume  isopleths  and  tower  wind  vectors 
from  the  Mt.  Iron  Report  Phase  I  Vol.  I. 
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Fig.  19   LINCOM  Case  90  (seabreeze  northwesterly)  wind  vectors 
shown  at  every  fourth  grid  point  (2km) .   Domain  is  25  x  40  km. 
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Fig.  20   LINCOM/RIMPUFF  and  Mt .  Iron  dose  isopleths  for  Case  90  at 
100m  resolution.   Domain  is  8.4  x  11.9  km. 
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Fig.  21   Mt.  Iron  Case  91  plume  isopleths  and  tower  wind  vectors 
from  the  Mt.  Iron  Report  Phase  I  Vol.  I. 
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Fig.     22       LINCOM    Case    91     (evening    northerly)    wind    vectors    shown    at 
every    fourth    grid    point    (2km) .       Domain    is    25    x    40    km. 
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Fig.  23   LINCOM/RIMPUFF  and  Mt.  Iron  dose  isopleths  for  Case  91  at 
100m  resolution.   Domain  is  8.4  x  11.9  km. 
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Fig.  24   Mt.  Iron  Case  110  plume  isopleths  and  tower  wind  vectors 
from  the  Mt.  Iron  Report  Phase  I  Vol.  I. 
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Fig.     25       LINCOM    Case    110     (Seabreeze    weak    westerly)     wind    vectors 
shown    at    every    fourth   grid    point    (2km) .       Domain    is    25    x    40    km. 
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LIN  COM  vs.    Mt.    Iron   Case    110 
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Fig.  26   LINCOM/RIMPUFF  and  Mt .  Iron  dose  isopleths  for  Case  110  at 
100m  resolution.   Domain  is  8.4  x  11.9  km. 
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Fig.  27   Scatter  plot  of  normalized  centerline  dose  exposures  for 
Mt.  Iron  versus  LINCOM/RIMPUFF.   Units  are  sec/nr 
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5.2   Discussion  of  Results 

As  shown  in  figs.  3-18,  the  LINCOM  generated  wind  fields  conform 
quite  closely  to  those  observed.  This  level  of  conformity  is  not 
usually  seen  in  diagnostic  models  and  is  mainly  due  to  the  "tower 
anchoring"  feature  discussed  in  section  3.1.  Since  the  towers  were 
not  located  at  exact  grid  points  and  tower  influence  decays  with  a 
modified  inverse  distance  squared,  the  more  exaggerated  flow 
divergence  indicated  by  the  towers  is  mitigated  on  the  gridded  wind 
field  by  LINCOM  dynamics.  However,  in  cases  90  and  91  this 
conformity  becomes  overzealous,  since  Tower  012  at  the  mouth  of 
Honda  Canyon  was  probably  in  error  (figs.  13,  15). 

As  an  example  of  the  "aliasing"  problem  discussed  in  section  3.1, 
LINCOM  is  also  overly  influenced  by  towers  at  the  bottoms  of 
canyons  where  the  winds  are  very  locally  controlled  but  do  not 
conform  to  the  bulk  of  the  boundary  layer  flow.  The  imposition  of 
strict  mass  conservation  as  in  LINCOM  1.1  should  retain  the  main 
flow  features  but  would  not  resolve  the  local  flows.  In  fact, 
tower  observations  will  generally  suggest  more  flow  divergence 
than  is  displayed  in  modeled  wind  fields  because  the  tower  winds 
respond  to  all  forcings,  including  the  ones  that  are  of  smaller 
scale  than  the  model  grid  resolution  such  as  turbulence  and  very 
local  terrain. 

This  leads  to  a  partial  dilemma  for  tracer  studies  in  hilly, 
complex  terrain.  As  at  Vandenberg,  available  roads  for  fixed  or 
mobile  samplers  are  often  along  the  bottoms  of  canyons.  If  the 
canyon  transects  the  plume,  canyon  samples  often  will  show  wider 
plume  widths,  reflecting  the  tendency  for  local  along-canyon  flow 
(see  LVDE,  1991) .  However,  most  of  the  plume  will  be  relatively 
unaffected  and  float  over  such  canyons.  Without  very  high 
resolution  and  complete  prognostic  dynamics,  numerical  models  may 
not  agree  well  with  data  when  the  sampling  does  not  conform  to  the 
bulk  of  the  flow.  Of  course  aircraft  data  would  tend  to  mitigate 
this  sampling  flaw. 

The  height  of  the  inversion  base  is  another  important  influence  on 
the  flow  field.  LINCOM/RIMPUFF  assumes  that  the  inversion  base 
parallels  the  terrain.  For  daytime  unstable  conditions,  this  is  a 
reasonable  assumption  as  shown  in  Kamada  et  al  (1990).  However, 
local  inversion  height  changes  are  suggested  in  cases  48  and  90, 
(see  section  2.4).  A  lowered  inversion  may  be  responsible  for  the 
high  winds  seen  at  the  Boathouse  in  case  90  (fig.  13).  The 
inversion  base  may  also  be  somewhat  flatter  than  terrain  parallel 
under  neutral  or  mildly  stable  conditions.  These  features  require 
a  complete  physics  package  and  cannot  be  resolved  by  the  dynamics 
of  simpler  diagnostic  models.  On  the  other  hand,  in  this  case  the 
"tower  anchoring"  feature  in  LINCOM  assists  modeling  accuracy  in 
the  winds  south  of  Honda  Ridge,  in  spite  of  the  crude  inversion 
height  assumption. 
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Figs.  19  -  26  and  the  corresponding  DICA,  FB ,  and  NMSE  scores 
suggest  that  LINCOM/RIMPUFF  matches  the  Mt.  Iron  isopleths  fairly 
well,  except  for  case  55  where  the  DICA  score  falls  below  0.2. 
However,  the  fractional  biases  indicate  that  the  LINCOM/RIMPUFF 
plumes  are  longer  than  Mt.  Iron  in  seven  of  the  eight  cases.   I.e., 
doses  are  higher  than  those  measured  at  longer  downwind  ranges, 
particularly  cases  28,  90,  and  110,  where  FB  are  less  than  -0.4, 
even  though  downwind  ranges  for  the  higher  dose  isopleths  are  quite 
comparable.    This  is  in  keeping  with  the  AFTOX/Mt.   Iron  and 
WADOCT/Mt.  Iron  comparisons  of  plume  centerline  dosages  (Kunkel  and 
Izumi,  1990;  Kunkel,  1991).   It  is  also  consistent  with  the  LVDE 
results. 

As  in  the  above  reports,  we  suggest  that  the  zinc  sulfide  wet 
aerosol  tracer  was  not  physically  inert  as  previously  assumed,  but 
probably  was  gradually  adsorbed  into  the  ground  canopy  as  the  plume 
was  advected  downwind.  This  is  consistent  with  oral  reports  that  by 
the  conclusion  of  the  Mt .  Iron  releases,  pools  of  dried  tracer 
material  surrounded  the  release  sites. 

For  most  cases,  the  DICA  and  NMSE  scores  yield  parallel  results. 
For  instance,  case  87  shows  the  best  fit  with  both  methods. 
However,  the  DICA  and  NMSE  scores  diverge  for  cases  31,  90,  and 
particularly  91.   The  NMSE  was  much  higher  for  case  91  than  the 
other  cases,  even  though  the  DICA  and  FB  scores  were  fairly  good. 
As  in  31  and  90,  the  NMSE  in  case  91  involves  a  large  number  of 
points  (column  2  of  Table  5.1.1).   However,  both  the  measured  and 
simulated  average  dosages  are  low  (columns  3  and  4),  as  is  the 
total  error  (column  5) .  But  the  NMSE  is  large  here  because  the  dose 
average  is  so  small.   Case  90,  with  a  much  larger  average  dosage 
and  error,  has  a  much  smaller  NMSE.  This  suggests  that  normalizing 
the  error  by  the  average  isopleth  value  may  skew  the  results  to 
favor  cases  with  large  average  dosages,  despite  having  large  errors 
prior  to  normalization. 

On  the  other  hand,  the  DICA  algorithm  does  not  seem  to  suffer  from 
this  problem.  It  sorts  the  LINCOM/RIMPUFF  model  results  for  the 
eight  cases  in  much  the  same  order  as  derived  by  a  qualitative 
visual  inspection,  but  more  accurately.  For  example,  a  visual 
inspection  (Fig.  23)  suggests  that  case  87  shows  the  best  fit.  It 
also  has  the  highest  DICA  score;  case  55  (fig.  22)  shows  the  worst 
fit  and  lowest  DICA  score.  However,  it  is  difficult  to  distinguish 
case  87  from  110  or  to  compare  two  rather  different  cases,  such  as 
28  and  91,  for  congruence  by  eye,  while  DICA  does  so  easily.  Given 
these  strengths,  we  suggest  that  DICA  be  utilized  as  a  standard 
performance  measure  where  applicable. 

Save  for  cases  55  and  90  (figs.  22  and  24),  the  modeled  plume 
thicknesses  are  comparable  to  the  Mt .  Iron  estimates.  In  cases  55 
and  90  RIMPUFF  overestimates  the  lateral  dose  footprint.  55  is  a 
special  case  discussed  in  more  detail  below.  However,  we  suspect 
that  the  problem  in  case  90  stems  from  underestimating  the  vertical 
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diffusion.  Wind/temperature  vertical  profiles  for  case  90  suggest 
near  neutral  conditions  with  significant  mechanical  turbulence. 
For  such  conditions,  RIMPUFF's  vertical  turbulence  parameterization 
will  tend  to  underestimate  the  extent  of  the  shear  based 
turbulence.  In  general,  it  is  also  hard  to  assure  accuracy  for 
tower  based  temperature  profiles  on  a  routine  operational  basis, 
due  to  the  need  for  freguent,  precise  sensor  calibrations.  Perhaps 
the  surface  radiation  (rather  than  profile  based)  stability 
algorithm  in  LINCOM  1.1  should  be  used  in  RIMPUFF  as  well. 

Whereas  the  releases  were  nearly  point  sources,  in  each  case  the 
modeled  plume  head  appears  thicker,  due  to  the  initial  100  meter 
puff  diameter.  A  smaller  initial  puff  diameter  could  have  been 
used.  So  this  is  not  a  wholly  necessary  model  artifact.  Part  of 
the  extra  thickness  is  also  due  to  the  interpolation  routine  used 
in  the  plotting  program,  as  seen  from  the  increased  thickness  in 
interpolating  the  hand  drawn  Mt.  Iron  plots.  However,  the  higher 
valued  dose  isopleths  are  generally  broader  in  RIMPUFF  than  the 
hand  drawn  isopleths.  This  is  not  as  true  of  the  lower  valued  dose 
isopleths,  since  diffusion  decays  exponentially  with  distance. 

We  also  see  that  whenever  the  plume  passes  over  Honda  Canyon, 
substantial  deviations  between  the  modeled  and  hand  drawn  isopleths 
can  result.   The  modeled  winds  are  influenced  by  towers  located  at 
the  bottom  of  Honda  Canyon,  while  in  cases  31,  48,  90,  and  91,  much 
of  the  real  plume  seems  to  have  floated  above  and  beyond  the 
sampling  lines  at  the  bottom  of  the  canyon.   The  LINCOM/  RIMPUFF 
tandem  cannot  reproduce  the  local  minimum  dosage  values 
seen  in  these  cases  at  the  bottom  of  Honda  Canyon.   Kunkel  (1990) 
reports  a  similar  modeling  constraint  for  the  WADOCT  model. 
Without   non-hydrostatic   capability,   i.e.,   both   buoyancy   and 
acceleration  in  the  vertical  momentum  eguation,  neutral  LINCOM' s 
ability  to  simulate  the  vertical  divergence  indicated  by  these 
cases  is  guite  limited.   For  example,  the  hydrostatic  version  of 
HOTMAC/RAPTAD  also  did  not  reproduce  the  local  minimums  seen  in 
cases  90  and  91,  even  though  HOTMAC  does  include  buoyancy  in  its 
vertical  momentum  eguation  (Yamada  and  Bunker,  1991) .   It  may  be 
that  resolutions  higher  than  500m  are  also  reguired  to  show  such 
effects. 

On  the  other  hand,  cases  31,  48,  and  90  (figs.  20,  21,  and  24) 
probably  indicate  aliasing  in  LINCOM  1.0.   Presumably,  the  mass 
conservation  integral  to  LINCOM  1.1  and  the  inclusion  of  buoyancy 
in  LINCOM  2.0  should  improve  the  simulation  in  such  cases. 

Perhaps  roughly  80%  of  the  accuracy  in  plume  dispersion  estimates 
depends  on  accurate  initial  wind  directions.  For  example,  case  28 
shows  that  a  few  degrees  error  (or  perhaps  real  wind  shift)  in 
specifying  wind  direction  can  lead  to  substantially  lower  DICA  and 
higher  NMSE  values. 

The  effect  of  downwind  towers  cannot  be  ascertained  for  the  weak 
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westerly  and  southwesterly  Seabreeze  cases:  28,  87,  and  110  (figs. 
19,  23,  and  26),  due  to  the  lack  of  data  beyond  the  sampling  line 
at  Arguello  Rd .  and  consequent  truncation  of  the  plumes. 

For  case  55  (fig.  22)  the  LINCOM/RIMPUFF  plume  is  much  longer  than 
the  measured  one.  As  discussed  in  the  case  flow  analysis  of 
section  2.3,  a  convective  cumulus  condition  with  vigorous  updrafts 
is  indicated  by  the  unstable  lapse  rates,  high  mid-April  inversion, 
and  patchy  cloud  cover  over  neighboring  sites  in  Honda  Canyon  and 
Scout  D.  Thus,  the  plume  may  simply  have  lofted  and  dispersed  over 
unsampled  areas.  Or  a  sustained  local  downdraft  may  have  forced 
adsorption  into  the  canopy.  The  former  interpretation  is  probably 
more  likely. 

There  are  data  problems  which  render  the  scatter  diagram  in  fig.  27 
less  meaningful  than  is  immediately  apparent.  1)  Since  we  compare 
only  eight  cases,  there  are  only  twenty-two  valid  data  points.  2) 
Only  the  hand-drawn  isopleths  were  available,  rather  than  actual 
bag  sample  exposure  levels.  Hence,  non-linear  interpolation  was 
needed  to  obtain  actual  values.  3)  The  sampling  lines  for  Mt.  Iron 
were  along  accessible  roads  rather  than  concentric  circles.  Thus, 
it  was  sometimes  difficult  or  impossible  to  compare  modeled  and 
measured  values  when  the  centerlines  diverged  in  direction.  4) 
Figure  19  seems  to  coincide  with  the  fractional  bias  results  which 
show  that  the  RIMPUFF  predictions  are  higher  than  those  observed  at 
lower  dosage  levels;  actually,  the  predicted  centerline  dosages  are 
higher  mainly  only  in  Honda  Canyon. 


6.   SUMMARY  AND  CONCLUSIONS 

Results  from  the  LINCOM/RIMPUFF  dispersion  model  were  compared  with 
eight  representative  cases  from  the  Mt.  Iron  tracer  study  conducted 
at  South  Vandenberg  AFB  during  1966.   We  conclude  that: 

1)  The  LINCOM  wind  fields  conform  fairly  closely  to  observations 
due  to  "tower  anchoring".  Tower  anchoring  allows  LINCOM  to  retain 
complex  flow  features  which  would  otherwise  require  the  use  of  a 
higher  resolution,  non-hydrostatic  prognostic  flow  model  with  a 
complete  physics  package.  However,  this  feature  is  a  mixed 
blessing  when  local  towers  misrepresent  the  main  flow.  For  tracer 
studies,  this  points  to  the  need  for  aircraft  sampling  in  addition 
to  ground  sampling  to  determine  local  as  well  as  bulk  flow 
features. 

2)  Plume  simulations  based  on  LINCOM/RIMPUFF  compare  fairly  well 
with  Mt.  Iron,  but  the  modeled  plumes  are  generally  somewhat 
longer,  as  is  consistent  with  other  published  reports.  Perhaps 
part  of  this  is  due  to  some  adsorption  of  the  zinc  sulfide  wet 
aerosol  onto  the  vegetation  canopy  during  transit. 

3)  Four  indicaters  were  used  to  compare  the  modeled  dispersion 
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pattern  with  observations.  The  most  useful  appear  to  be  scores 
from  the  dose  isopleth  correlation  area  (DICA)  method  and  the 
fractional  bias  (FB)  .  We  also  used  the  normalized  mean  square 
error  (NMSE)  and  a  centerline  scatterplot  comparison.  The  NMSE 
seems  to  suffer  from  a  tendency  to  overweight  errors  when  dosages 
are  relatively  low.  The  centerline  scatterplot  involves  too  few 
data  points  and  suffers  the  effects  of  some  other  Mt.  Iron  data 
idiosyncrasies  discussed  in  section  5.2. 

4)  Lateral  plume  thicknesses  compare  fairly  well  with  those 
observed,  except  perhaps  when  turbulence  intensities  were  mis- 
calculated, due  to  observational  error,  limited  spatial  resolution, 
or  deficiencies  in  the  turbulence  parameterization  algorithms. 

5)  The  low  dosages  recorded  in  Honda  Canyon  indicate  that  plumes 
tend  to  float  over  the  canyon  rather  than  parallel  the  terrain 
contours.  LINCOM/RIMPUFF' s  ability  to  handle  this  feature  is 
rather  limited,  since  a  complete  vertical  momentum  equation  is  not 
included.  Other  diagnostic  and  prognostic  model  simulations  of  Mt. 
Iron  have  also  shown  an  inability  to  reproduce  such  local  dosage 
minima.  Resolutions  higher  than  500m  may  be  helpful  in  refining 
such  simulations. 
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8.  APPENDICES 

8.1    LINCOM  Model  Theory 

8.1.1  Introduction 

We  simulated  plume  dispersion  for  each  of  the  above  ten  flow  types 
with  a  combined  flow/puff  model.  LINCOM,  the  flow  model,  extends 
linear  neutral  potential  flow  theory  (over  a  single  hill)  to  a  more 
general  case  involving  mesoscale  complex  terrain.  LINCOM  is  a 
diagnostic  code  which  employs  highly  efficient  solutions,  which  are 
spectrally  based  in  the  horizontal  plane  and  analytic  in  the 
vertical.  A  "look-up  table"  of  pre-computed  results  streamlines 
this  procedure  even  further.  Based  on  data  from  ten  or  more 
towers,  LINCOM  1.0  can  output  a  100  x  140  gridded  wind  field  for 
five  vertical  levels  in  the  boundary  layer  at  500m  horizontal 
resolution  in  less  than  one  minute  on  a  386/33  PC. 

The  flow  model  is  initialized  with  mean  uniform  velocity,  pressure, 
and  potential  temperature  fields,  (U0,  V0,  W0,  and  P0)  ,  using  data 
from  each  tower  one  at  a  time.  Thus,  at  each  tower  the 
perturbation  quantities  u,  v,  w,  and  p  induced  by  the  terrain  are 
all  zero.  LINCOM  then  solves  the  governing  hydrodynamic  equations 
solely  in  terms  of  the  perturbation  quantities;  and  adds  the 
perturbation  fields  to  the  mean  fields  to  obtain  the  final  results: 
U  =  U0  +u,  V  =  V0  +  v,  etc.  So  at  any  distance  away  from  that  tower 
the  u,  v,  w,  and  p  can  assume  non-zero  values.  We  repeat  this 
procedure  for  all  n  towers  to  obtain  n  separate  initial 
perturbation  wind  fields.  These  fields  are  combined  with  a  modified 
1/r2  weighting  scheme  which  warps  each  tower's  circle  of  influence 
into  ovals  that  are  essentially  congruent  with  the  elevation 
contours  for  the  surrounding  terrain.  This  method  ensures  that  the 
final  field  will  exactly  match  the  tower  data,  while  tower 
influences  extend  strongly  only  along  similar  terrain  elevations. 
In  this  way,  even  under  differing  stability  conditions,  a 
sufficient  density  of  towers  will  tend  to  restore  some  of  the 
non-neutral  physics  and  force  the  final  field  to  closely  resemble 
the  measured  flow. 

One  problem  with  this  procedure  is  that  mass  continuity  can  be 
compromised.   Thus,  in  LINCOM  1 . 1  we  instead  find  the  set  of  mean 
fields  which  provides  the  best  fit  to  the  existing  set  of  tower 
data  when  added  to  the  resultant  perturbation  fields. 

8.1.2  Governing  equations 

To  discuss  the  modeling  ideas  in  more  detail,  we  begin  with  a 
truncation  of  the  full  atmospheric  equation  set  which  neglects 
density  variations  due  to  effects  other  than  gravity.    This 
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standard  Boussinesq  set  forms  the  basis  for  almost  all  atmospheric 
modeling.  Included  are  the  three  dimensional  momentum  equations 
(neglecting  horizontal  shear) : 

du     du     du     du  1  Up         1  d 

—  +  U —  +  V —  +  W fv= t  ,    (8.1.1a) 

dt     dd     dy     dz  p    dx        p    dz  xz 

dv    dv    dv    dv         1  dp   1  d 

—  +  U —  +  V —  +  W —  +fu  = t      ,    (8.1.1b) 

dt  dx  dy     dz  p    dy        p    dz   yz 

dw     dw     dw     dw      1  dp    1  d  g<59 

—  +  U —  +  V —  +  W —  = t   -   ,    (8.1.1c) 

dt  dx  dy  dz  p    dz         p    dz    zz      6 

the  thermodynamic  energy  (or  potential  temperature)  equation, 

de         de         de    de  l  de 

—  +  U —  +  V —  +  W —  +  Wr  = ;  (8.1.2) 

dt  dx  dy     dz  p    c  dz 

P 

the  shallow  convection  form  of  the  mass  conservation  or 
incompressible  fluid  continuity  equation  (applicable  to  boundary 
layer  models) ,  and  the  ideal  gas  equation  of  state, 

du    dw        dw 

—  +  —  +  —  =0      ;      p  =  pRT      .  (8.1.3) 
dx         <3y    dz 

Here,  f,  the  Coriolis  parameter  is  defined  as  2flsin  e,  where  n  is 
the  earth's  rate  of  rotation.  r,  the  potential  temperature  lapse 
rate,  is  dQ/dz.  R  and  c  are  the  gas  and  specific  heat  constants 
for  air,  and  the  t-  represent  eddy  stresses.  We  wish  to  simplify 
these  equations  because  their  complete  solution  at  all  relevant 
scales  of  motion  in  what  is  termed  a  direct  simulation  is  beyond 
the  forseeable  capacity  of  computers.  That  is,  the  current  direct 
simulation  limit  on  CRAY  class  supercomputers  is  "  107  grid  points 
for  mesoscale  flows.  A  3-D  direct  simulation  would  require  "  1023 
grid  points  to  resolve  the  larger  mesoscale  (105m)  circulations, 
the  dissipation  of  turbulent  kinetic  energy  which  occurs  down  to 
10"3m,  and  the  intervening  scales  of  turbulent  and  forced  motions. 
Since  operational  models  use  less  than  105  grid  points,  sub-grid 
scale  effects  must  be  parameterized  rather  than  modeled  directly. 
Much  of  the  energy  is  resolved  on  the  grid  itself;  however, 
sub-grid  motions  cannot  be  neglected  because  the  non-linear 
eddy  stress  dynamically  links  all  scales  of  turbulence  in  an  energy 
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cascade.   Thus,  without  dissipation  into  internal  energy,  that  is, 
heat,  the  modeled  turbulence  kinetic  energy  would  simply  grow 
without  bound. 

Other  difficulties  abound.   For  example,  due  to  the  continuum  of 
turbulence  scales,  simple  analytic  forms  for  eddy  stress  do  not 
exist.   The  better  approximations  are  spatially  dependent.   They 
also  involve  higher  order  correlations  between  both  supra  and 
sub-grid  velocity  and  pressure  fluctuations  whose  magnitudes  are 
still  being  debated.   Non-linear  advection  also  creates  annoying 
cross  terms  which  are  usually  simply  ignored  when  the  equations  are 
written  in  non-orthogonal  terrain  following  coordinates.   Thus, 
even  with  supercomputers,  the  full  or  "primitive"  equation 
mesoscale  models  are  still  too  unwieldy  for  operational  purposes. 

Thus,  to  simplify  in  a  way  which  hopefully  blends  suitability  with 
convenience,  we  can  make  linearizing  approximations  of  the  sort, 

UdU/dx  =  (U0  +  u)d(U0  +  u)/3x  "  U0du/3x    .  (8.1.4) 

Fourier  transforming  the  linearized  equations  to  spectral  space 
will  then  allow  symbolic  matrix  inversion  solutions  for  each 
wavenumber  for  the  entire  domain  at  once,  rather  than  grid  point  by 
grid  point  as  in  finite  differencing  methods.  However,  we  must 
remain  aware  that  the  linearization  procedure  neglects  the  higher 
order  products  of  perturbed  quantities,  such  as  udu/dx.  This  is 
only  valid  if  udu/dx  <<  U0du/dx. 

To  suit  this  level  of  modeling  we  match  the  above  approximation 
with  a  first  order  approximation  of  the  eddy  stress,  i.e.,  ruz  =  -K 
du/dz,  where  K  is  the  horizontally  diffusivity  coefficient. 

Moreover,  for  stationary  flows  or  at  least  if  dV/dt  <<  UdU/dx,  we 
can  also  neglect  time  derivatives.  This  inequality  implies  that 
1/T  <<  u/L,  where  T  and  L  are  characteristic  time  and  length  scales 
for  flow  changes.  For  a  typical  T  "  6hr  and  u  "  lm/s,  we  need  an 
L  <<  20  km.  For  u  "  lOm/s,  we  need  an  L  <<  200  km.  These  are  both 
reasonable  for  Vandenberg. 

In  LINCOM  1.0  and  1 . 1  we  also  neglect  buoyancy  forces  caused  by  an 
uneven  density/temperature  field.  So  neutral  LINCOM  cannot 
diagnose  thermally  induced  slope  flows  or  seabreezes,  except  as 
indicated  within  the  mean  flow  given  by  the  towers  or  the  later 
tower  based  objective  analysis.  this  neglect  allows  a 
decomposition  of  the  flow  field  into  orthogonal  components  which 
are  pre-calculable  for  a  given  domain  and  thus  lead  almost 
immediately  to  a  final  wind  field.  Also,  thermal  forcing  at  scales 
less  than  10  km  is  often  subordinate  to  the  dynamic  pressure  forces 
which  neutral  LINCOM  does  account  for. 
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At  any  rate,  we  are  left  with  the  following  simplified,  steady 
state,  linearized  equation  set  for  the  perturbed  part  of  the  flow: 


du 


dx 


du 


dy 


dp  d2u 

—  +  fv  +  K  

dx  dz2 


(8.1.5a) 


dv 
dx 


u0  —  +  v0 


dv 

ay 


aP       a2v 

fu  +  K  

dy  dz1 


(8.1.5b) 


dw 


dx 


dw 


ay 


dp 
dz 


+  K 


d\ 


dz- 


(8. 1.5c) 


The  continuity  and  ideal  gas  equations  remain  unchanged.  The  1/p 
density  factor  has  been  absorbed  into  the  perturbation  pressure,  p. 
As  discussed,  we  have  dropped  the  buoyancy  term,  g<59/9.  We  have 
also  dropped  the  small  vertical  terms.  The  solution  set  is 
analytic  for  any  choice  of  the  perturbation  quantities.  (u,  v,  w, 
p,  and  K) 


8.1.3 


Solution  method 


To  solve  eqns.  (8.1.2,  3,  and  5)  we  drop  the  Coriolis  term  and 
begin  with  a  simplified  2-D  version  of  the  model  in  order  to  sketch 
the  main  ideas.  Later,  we  add  more  complex  refinements  included  in 
the  real  3-D  model.  First,  we  Fourier  transform  the  x  component  of 
velocity  into  wave  number  (k)  space  in  the  fashion, 


no 


u(k,z)  =  \        u(x,z)  exp(ikx)  dx 


(8.1.6a) 


u(x,z) 


27T 


00 


j   u(k,z)  exp(-ikx)  3k 

—  00 


(8.1.6b) 


where  eqn.  (8.1.6b)  is  the  back  transform  of  u(k,z).   Note  that  the 
z  direction  does  not  participate  in  the  transform  and  is  written 
only  to  remind  us  of  the  height  dependence  of  u.   Then  by 
differentiation  of  eqn.  (8.1.6a), 
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du        w 

—  =  ikf   u(k,z)  exp(-ikx)  dk     =  iku    .  (8.1.7a) 

dx       -oo 


Thus,  the  momentum  equation  for  u  may  be  written  as, 

00  00 

U0ikj   u(k, z) exp(-ikx) 3k  +  ikj  p(k, z) exp(-ikx) 3k  - 

—00  —00 

3    oo 
K  —  \   u(k, z)exp(-ikx) 3k  =  0.      (8.1.7b) 

dz  -oo 

So  the  sum  of  the  wave  number  components  must  total  to  zero.  For 
each  wave  number,  the  2-D  version  of  the  governing  equations  will 
then  look  like, 

ikU0u(k,z)  +  ikp(k,z)    -  K  32u/3z2  =  0      ,  (8.1.8a) 

ikU0w(k,z)  +  3p(k,z)/3z  -  K  32w/3z2  =  0     ,         (8.1.8b) 

iku(k,z)    +  3w(k,z)/3z  =  0     .         (8.1.8c) 

By  virtue  of  the  Fourier  transform,  the  sum  of  the  perturbation 
quantities  at  each  wave  number,  k,  summed  over  all  k,  is  equivalent 
to  their  determination  for  each  grid  point  along  x  in  real  space. 
However,  it  remains  to  determine  the  dependence  of  these 
perturbation  quantities  along  the  z  axis.  Thus,  eliminating  w(k,z) 
and  p(k,z)  from  this  set  leads  to  a  fourth  order  ordinary 
differential  equation  (O.D.E.)  with  constant  coefficients, 

d4u(k,z)  32u(k,z)     ikU0k2u(k,z) 

-   (k2  +  ikU0/K) +  =  0. 

3z4  dz2  K 

(8.1.9) 

As  is  the  usual  fashion  for  ordinary  differential  equations,  we  can 
assume  exponential  solutions  in  z  of  the  form,  exp(-orz),  where  a   is 
a  complex  number.   This  leads  to  the  fourth  order  polynomial 
equation, 
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a4    -    (k2   +    ikU0/K)cr   +    ikU0k2/K   =    0 


(8.1.10) 


which  is  characteristic  of  the  fourth  order  O.D.E.   This  equation 
has  the  four  roots, 


a   = 


+/-    |k|  ,       and 

+  /-  \/(ikU0/K) 


(8.1.11) 


However,  the  two  positive  roots  lead  to  unlimited  increases  in  u 
with  z.   Thus,  only  the  two  negative  roots  represent  components  of 
the  real  solution, 


u(k,z)   =   U,  exp(v/(ikU0/K)  z)  +  U2  exp(-|k|z) 


(8. 1. 12a) 


or 


u(k,z)   =   U,  exp(-a,  z)  +  U2  exp(-a2  z) 


(8.1. 12b) 


Note  for  each  component  solution  that 


du(k,z)/3z  =  — or-U-  exp(-cr  z) 


(8.1.  13) 


So  in  general  a  is  equivalent  to  the  differential  operator, 


a   =  d/dz 


(8.1.14) 


We  will  use  this  below  in  determining  U,  and  U2.  But  before  doing 
so,  we  focus  on  the  exponential  terms.   First,  we  introduce  the 
inner  and  outer  length  scales,  £,  and  L.   In  the  simplest  case,  L 
is  the  horizontal  dimension  for  a  solitary  hill.    This  will 
introduce  only  one  wave  number,  k  =  1/L,  into  the  flow  response. 
We  define  the  inner  scale,  £,  by 


£  In  (£/z0)  =   kzL 


2T  ^ 


(8.1.15) 


where 


is   the  roughness   length   (Jackson  and   Hunt,   1979) 


Physically,  we  are  assuming  that  the  hill  interacts  strongly  with 
the  region  below  £  through  turbulent  exchange.  Thus,  below  £  is  a 
region  of  large  wind  shear,  whereas  above  £  turbulent  exchange  is 
small  and  resulting  in  nearly  inviscid  potential  flow. 
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For  neutral  flow  Monin-Obukhov  similarity  theory  specifies  for 
horizontally  homogeneous  surface  layers  that 


U0(z)  =  u»ln (z/z0)  /k  ,   and  (8.1.16a) 

K(z)   =  u*z//<   ,  (8.1.16b) 


where  k,      the  von  Karman  constant   is  =  0.4.    Thus,   U0/K  - 

'■o. 


In (z/z0)  / k2 .      This  leaves, 


a,  =v/(ikU0/K)  =  (l  +  i)v/(ln(z/z0)/k)  A/2  ,  (8.1.17) 

where  we  have  used  Vi  =  (l  +  i)/v/2.   Now  at  the  height  z  =  E,  we  have 

a,  =  (l  +  i)/v/(2£)    ,  (8.1.18) 

while  the  second  root  is  simply, 

tt1  =  -Ikl  =  -l/L    .  (8.1.19) 


Thus,  we  have  framed  the  advection  velocity  component  at 
wavenumber,  k,  in  terms  of  the  two  Jackson-Hunt  length  scales. 

Returning  to  the  U,  and  U2  coefficients,  for  real  flows  with  a 
non-slip  bottom  boundary,  u(k,z)  must  become  zero  at  z  =  0.  From 
eqn.  (8.1.12b)  this  forces 

U,  =  -U2    .  (8.1.20) 

For  the  magnitude  of  U, ,  we  assume  to  first  order  that  the  vertical 
wind  near  the  surface  is  given  by  the  vertical  component  of  the 
wind  vector  parallel  to  the  surface,  whose  speed  is  taken  to  be 
equal  to  the  horizontal  wind  speed.  That  is,  U  •  Vh  =  w(0),  where 
h  is  terrain  height  and  U  is  the  steady  state  wind.  By  the  Fourier 
transformation,  h(x)  becomes  -ikh(k)  ,  where  h(k)  is  defined  by 


1   oo 

h(x)  =  j   hk  exp(-ikx)dk    .  (8.1.21) 

2rr    -oo 
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This    leads    to    the   condition, 

W,  w2 

-ikh(k)    =  +      .  (8.1.22) 

U01  U02 

Here  we  have  approximated  the  mean  wind  U0  from  eqn.  (8.1.16a)  in 
terms  of  two  components,  U01  and  U02,  as  functions  of  a,  and  a2  and, 
hence,  £  and  L.  Namely,  U0,  and  U0~,  are  the  mean  wind  velocities  at 
z  =  £  and  z  =  L.  We  approximate  K,  the  turbulent  diffusivity, 
similarly  in  terms  of  t  and  L.  This  implies  that  a  perturbation 
penetrating  to  larger  heights  will  be  subject  to  larger  mean 
velocities  and  dif f usivities  than  less  penetrative  perturbations. 

At  any  rate,  from  eqns.  (8.1.8)  and  (8.1.14),  we  have 

w(k,z)=iku(k,z)/a     .  (8.1.23) 

Combining    eqn.     (8.1.20)    with    (8.1.22)    and    (8.1.23)    gives 

ik  ik 

+  u,    =   -ikh(k)         .  (8.1.24) 


alU01  a2U02 

With     LU0?     >>     £U01     we     can     neglect     the     first     term     and     with     eqn. 
(8.1.19); 

U,    «    -|k|U02   h(k)  .  (8.1.25) 


Substituting  this  into  eqn.  (8.1.12),  we  finally  obtain  for  the 
simple  2-D  case  the  velocity  perturbation  for  a  single  wave  number, 
k,  at  height  z  above  the  terrain, 


h(k) 

u(k,z)  =  U0(L)[e-7/L  -  e(,+l)zA/2,)]   .  (8.1.26) 

L 
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We  then  back  Fourier  transform  to  real  space  and  add  the 
perturbation  to  the  mean  velocity  to  produce  the  total  U  field.  A 
plot  of  the  vertical  profile  shows  that  the  maximum  perturbation  in 
2-D  occurs  at  (3  /2V2)£  "  3.32.  That  is,  hills  induce  pressure 
gradients  which  interact  with  the  shear  stress  to  force  advection 
to  produce  a  low-level  jet,  as  in  Jackson-Hunt  theory. 

Generalizing  to  3-D  flow  requires  that  we  re-introduce  the  Coriolis 
term  and  add  an  equation  and  terms  for  the  y  component.   In  the  2-D 
case  we  solved  the  O.D.E.  system  by  eliminating  variables  and 
consolidating  equations  until  only  one  fourth  order  O.D.E. 
remained.  However,  for  3-D  cases  this  procedure  is  too  cumbersome. 
So  we  re-organize  the  mathematics  by  writing  the  governing 
equations  in  matrix  form  and  solve  them  systematically  through 
elementary  row  operations  on  the  matrix,  rather  than  by 
elimination. 

For  the  3-D  case  the  form  of  the  Fourier  transform  and  back 
transform  now  become, 


oo  oo 

u(k,m,z)    =      J         \        u (x, y , z) exp(ikx   +    imy) dx   dy      ,         (8.1.27a) 

—oo        -oo 


1        oo  oo 

u(x,y,z)    =  —   J        J      u(k,m, z) exp(-ikx-imy) 3k  dm      ,         (8.1.27b) 

27T     -oo        -oo 


or  more  precisely,  using  the  discrete  (fast)  Fourier  transform 
pair, 


n-1  n-1 
u(k,m,z)  =   £    £   u(x,y,z)Wnkx  Wnmy    ,  (8.1.28a) 

k=0  m=0 


1  n-1  n-1 

u(x,y,z)  =  2   £   u(k,m,z)Wn.kx  Wn.my    ,         (8.1.28b) 

n2  x=0  y=0 

where  W  =  exp(i27r/n),  and  similarly  for  v,  w,  p,  and  h.   Then  the 
momentum  and  continuity  equations  for  the  Fourier  transformed 
components  at  each  wave  number  pair,  (k,m) ,  may  be  written  in  the 
form  of  the  matrix  equation, 
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0 

w(k,m) 

0 

p(k,m) 

0 

(8. 1.29) 


where,  using  eqn.  (8.1.14), 


C  ■  i(kU0  +  mVQ)  -  Ka' 


(8.1.30) 


Rather  than  by  elimination,  we  can  solve  this  set  systematically 
through  elementary  row  operations  (cf.  any  text  on  O.D.E.  systems) 
In  place  of  eqn.  (8.1.10)  the  characteristic  equation  becomes, 


C(k2  +  a:) 


(8. 1.31) 


We  replace  eqn.  (8.1.15)  for  £  by 


£  In  (l/z0)  =  rL/cos  6 


(8.1. 32) 


where  the  mean  wind  vector,  V  =  (U0,V0)  ,  L  =  1/k,  and  6  is  the  slope 
angle,  such  that 


|V| | (k,m) |cos  6  =  V  •  (k,m) 


(8.  1.33) 


The  surface  boundary  condition  for  w  changes  from  eqn.  (8.1.22)  to 


-ikh 


k,m 


W, 


W- 


(k,m)  •  V01 


(k,m)  •  v02 


(8.1.34) 


For  winds  not  perpendicular  to  a  2-D  hill,  Troen  (1986)  showed  that 
eqn.  (8.1.26)  still  gives  the  perturbation  of  the  component  normal 
to  the  hill.  That  is,  the  first,  outer  scale,  term  in  L  is  not 
affected  by  wind  direction,  but  £,  given  by  eqn.  (8.1.32),  is  now 
larger  for  off-normal  winds;  below  £,  the  wind  speed  perturbation 
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is  no  longer  given  by  a  simple  cosine  response. 

An  idealized  neutrally  stable  boundary  layer  should  have  no  abrupt 
top.  In  reality  the  top,  H,  is  usually  specified  by  other 
conditions,  such  as  a  subsidence  inversion  (chronic  for  Vandenberg) 
or  the  advected  "fossil"  remains  of  a  capping  inversion  formed  when 
the  upstream  boundary  layer  was  last  convectively  unstable.  So  we 
must  account  for  H.  We  also  expect  that  LINCOM  will  be  used  in 
cases  where  the  atmosphere  is  not  truly  neutral.  That  is,  the 
inverse  Obukhov  length,  1/Lmn  will  not  be  zero.  In  that  case,  Lmo 
can  be  measured  from  surface  data  and  used  in  similarity  functions 
in  order  to  obtain  the  vertical  profiles  of  mean  wind  speed  and 
diffusivity.  E.g.,  the  mean  wind  and  diffusivity  profiles  may  be 
obtained  from  the  following  extended  boundary  layer  similarity 
profiles  of  Sorbjan  (1989).   He  suggested  that 


V0  =  2.5u*(ln(z/z0)  +  5  z/L  -  z/H  -  2.5  z~/HLmo)  )   ,  (8.1.35a) 

K  =  0.4  U*Z(1  -  z/H)/(l  +  5  z/Lmo)  ,                 (8.1.35b) 
for  stable  cases  and 

V0  =  2.5u,(ln(z/z0)  +  <pm       ,  (8.1.36a) 

K  =  0.4  u,z(l  -  28  z/Lmo)1/4   ,  (8.1.36b) 
for  unstable  cases.   Here  we  supply  a  new,  more  efficient  algorithm 

0m  =  (1.19  +  0.23*ln(-z/Lmo)  )2   .  (8.1.37) 


8.1.4   Obukhov  length 

Unless  sensors  are  frequently  calibrated,  potential  temperature 
differences  between  vertical  levels  cannot  be  measured  accurately 
enough  on  a  routine  operational  basis  to  maintain  confidence  in 
atmospheric  surface  layer  stability  estimates.  Thus,  LINCOM  1.1 
does  not  determine  Obukhov  length  from  the  potential  temperature 
and  wind  speeds  differences  between  vertical  levels.  Instead,  it 
relies  on  measured  solar  and  thermal  radiation  values  or  estimates 
thereof.   This  type  of  method  is  useful  for  Vandenberg  because  the 
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near  coastal  areas  are  chronically  shrouded  with  stratus  cover, 
while  the  sky  is  often  quite  clear  several  kilometers  inland.  Thus, 
for  diagnostic  models  such  as  LINCOM,  Obukhov  lengths  need  be 
computed  for  only  a  few  areas  between  sun  and  shade,  based  on 
moderate  differences  in  roughness  length  and  local  wind  speeds. 

We  also  included  a  routine  to  estimate  downwelling  solar  and 
thermal  radiation,  if  measured  data  is  not  available.  This  routine 
requires  a  fractional  cloud  coverage  estimate  as  well  as  other 
commonly  available  input  data,  such  as  inversion  height,  screen 
level  wind  speed,  temperature,  and  relative  humidity. 

All  solar  radiation  models  require  current  date  and  time  to  compute 
sun  angle,  \p .  This  is  done  in  LINCOM  using  standard  algorithms.  We 
modified  a  simple  ASHRAE  model  for  downwelling  solar  direct  and 
diffuse  radiation  (Iqbal,  1983)  by  computing  sine  curve  fits  to  the 
seasonal  adjustments  for  apparent  extra-terrestrial  radiation, 
atmospheric  optical  air  mass,  and  column  length  of  precipitable 
water,  W. 

If  cloud  base  height,  zc,  and  boundary  layer  height,  H,  are  not 
available,  the  solar  transmissivity  through  clouds,  tc,  is 
estimated  using  the  simple  algorithm, 


tc  =  (1  -  0.6FCC)   ,  (8.  1.38) 

where  FCC  is  fractional  cloud  coverage  of  the  sky.  FCC  can  be 
obtained  from  ground  based  observation  and/or  on-line  GOES 
satellite  data,  using  the  algorithm  described  in  Skupniewicz  et  al, 
(1991).  If  zc  and  H  are  available  (from  rawinsonde  data  or 
aircraft  landing  reports) ,  the  solar  transmissivity  for  a  non- 
reflective  ground  surface  can  be  estimated  more  accurately  from  the 
method  of  Liou  and  Wittman  (1979).  This  method  uses  sun  angle  and 
column  height  of  precipitable  water  within  the  cloud  in  a  bivariate 
polynomial  regression  based  on  results  from  an  accurate  multi- 
stream  discrete  ordinates  model.  Currently,  we  assume  that  the 
cloud  coverage  is  stratiform  and  confined  to  the  boundary  layer, 
with  a  liquid  water  content  of  0.78  grams/meter  of  cloud  depth. 
Hence,  the  only  inputs  required  are  date,  time,  cloud  base  height, 
and  boundary  layer  depth..  The  algorithm  can  be  extended  easily  to 
include  other  cloud  types  and  water  content.  The  regression  form  is 

3      3 

Tc(fi0,W)  =  S    E  h,    /i'o  W>    ,  (8.1.39) 

i=0   j=0  J 


where  iiQ   is  solar  zenith  angle,  W  is  precipitable  water,  and  the  b- 
are  the  coefficients  obtained  from  the  regression.   For  actual  non- 
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zero  surface  albedoes,  As  (default  value,  0.15),  we  modify  the  cloud 
solar  transmissivity ,  using  the  algorithm  of  Kamada  (1984), 


i 


FCC 


1  -  A  A 

rcm  =   — (1  -  d0)Tc  +  d0) 

(1  -  0.12AC)  (1  -  0.1AS) 
L  J 

(8.1.40) 

where  Ac  is  cloud  top  albedo  (  default  value  for  stratiform  clouds 
is  0.55),  d  =  0.001068,  and  <p  is  in  degrees  latitude.  The  total 
downwelling  solar  radiation  is  then 

SOLi    =    (I0sin(^)     +    Idirf)Tcm         ,  (8.1.41) 

where  Idiff  is  the  downward  diffuse  solar  radiation  at  the  surface. 

The  downwelling  thermal  radiation  computation  is  initiated,  using 
the  algorithm  of  Martin  and  Berdahl  (1983),  which  employs  surface 
dewpoint  temperature,  Td  ,  hour  of  the  day,  Hr,  and  pressure,  p,  to 
estimate  the  effective  clear  sky  emissivity, 

€cs  =  0.711  +  0.0056  Tdp  +  0.00073  Tdp2  + 

0.013cos(0.262Hr)  +  0.00012(p  -  1000)   .       (8.1.42) 

Td  is  readily  obtained  from  the  relative  humidity  and  temperature, 
using  standard  formulas.  The  emissivity  is  then  modified  for  clouds 
according  to  cloud  base  height,  zc,  and  fractional  cloud  coverage, 
FCC.   Thus,  we  have 

€c  =  ecs  +  0.85FCC(1  -  ecs)  exp  ( 1 .  22xl0"4  Zc)   .        (8.1.43) 

Again,  boundary  layer  stratus  clouds  are  assumed  here  but  other 
cloud  types  are  readily  included.  Downwelling  thermal  radiation  is 
computed  by  assuming  the  cloudy  or  clear  sky  to  be  a  grey  body 
thermal  emitter,  such  that 

IR1  =  eca94    ,  (8.1.44) 
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where  a     =     5.67x10   is  the  Planck  black  body  constant.  Total 

downwelling  radiation  at  the  earth's  surface  is  obtained  by 
combining  solar  and  thermal  contributions  via, 

RADi  =  (1  -  A) SOLI  +  IR1     .  (8.1.45) 


With  the  radiation  component  of  the  surface  energy  budget  computed, 
we  can  obtain  the  atmospheric  stability  and  temperature  flux  from 
the  ground  surface  to  the  air.  As  a  measure  of  atmospheric 
stability,  the  Obukhov  length  is  defined  as, 

L  =  u*2e/g*e«    ,  (8.  1.46) 

where  u.  is  the  surface  layer  friction  velocity,  e  =  T( 1000/p)0  285 
is  the  near  surface  potential  temperature,  g  is  gravitational 
acceleration,  k  is  0.4,  the  von  Karman  constant,  and  G*  is  the 
Obukhov  temperature  scale.  The  temperature  flux  can  be  defined  as 
the  statistical  correlation  between  vertical  velocity  and  potential 
temperature  perturbations,  and  is  also  given  by 


w'e'0  =  -u*6*    .  (8.1.47) 

Thus,  given  the  downwelling  radiation,  we  can  iterate  between 
estimates  of  L  and  estimates  of  surface  temperature  flux  until  both 
guantities  converge.  This  reguires  that  we  compute  both  u*  and  e*. 
u*  comes  from 


u.  =  U(z) 


Ln(z/z0)  -  tm  (8.1.48) 


where  U(z)  is  the  mean  windspeed  at  height  z,  and  z0  is  the  surface 
vegetative  canopy  roughness  length  (typically  —  1/7  the  mean 
vegetation  height).  <pm  is  given  by  egn.  (8.1.37).  Analogous  to  u,, 
6*  is  given  by 


9,  =  69(z)l k ,  (8.1.49) 

Ln(z/z0)  -  *h 
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From  Dyer  and  Bradley  (1982),  we  have 


%   =  2  Ln(  h     (1  -  J  1  -  14z/L  ) )    .  (8.1.50) 


To  obtain  the  ground  surface  "skin"  temperature,  e00,  we  assume  that 
the  potential  temperature  difference  between  the  roughness  height, 
z0,  and  height,  z,  is  given  by  68,  and  that  the  temperature 
difference  across  the  laminar  layer  between  z0  and  the  surface  is 
given  by  9,.   Thus, 

600  =  9  -  (9,/k)  (Ln(z/z0)  -  %)    -   9,    .  (8.1.51) 

Since  69  is  not  known,  9*  is  initially  set  to  zero  and  the  skin 
temperature  is  set  equal  to  the  screen  level  potential  temperature. 
The  temperatures  then  diverge  toward  equilibrium  values  by 
iteration.  The  net  radiative  budget  at  the  surface  is  given  by 

NETRAD  =  RAD1  -  eg  a  9004     ,  (8.1.52) 

where  the  ground  surface  emissivihy,  e  ,  has  a  default  value  of 
0.95.  Following  the  Penman-Monteith  model  (1948,  1965), 
temperature  flux  into  the  ground  is  estimated  as,  Q  =  0.15  NETRAD 
with  the  stipulation  that  for  stable  conditions  (if  L  >  0) ,  Q  is 
3.3  times  larger.  The  Penman-Monteith  equation  is  used  to  estimate 
the  temperature  flux, 


-(7  (-NETRAD  +  Q  )  -  w'q'0) 


W  '  9  '  0  = -  0  .  8  4  w  •  q  '  0 , 

(8. 1.53) 


(  QCP  (xgscc  +  7)  ) 


where  w'q'0  is  the  humidity  flux,  p  is  air  density,  c  =  1005  J  kg"1 
K"1  is  the  heat  capacity  of  air,  X  is  soil  relative  humidity,  scc  is 
the  change  rate  of  specific  humidity  with  temperature  for  saturated 
air,  and  7  =  cp/Lv  =  0.0004°K"1  is  the  psychrometric  constant.  In 
turn  these  latter  variables  are  obtained  from  standard  algorithms. 
The  Obukhov  temperature  scale  is  obtained  from, 

9*  =  -w'e'o/u*    .  (8. 1.54) 
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In  this  scheme,  note  that  under  neutral  conditions  when  the 
temperature  flux  falls  to  zero,  9*  will  vanish  and  L  becomes 
infinite.  Thus,  to  avoid  infinities  during  iterative  numerical 
evaluation,  it  is  better  to  compute  the  inverse  Obukhov  length, 
1/L. 

From  a  rough  bivariate  analysis  of  our  results,  we  found  that  a 
useful  first  guess  is 

1/L  =  0.000674(300  -  RADi)  Ln ( 10  Z /  z0)  /U2  ( z)  .     (8.1.55) 

In  order  to  cover  a  wide  range  of  radiation  values,  wind  speeds, 
temperatures,  and  roughness  lengths,  the  iteration  procedure  must 
be  guite  robust,  otherwise  convergence  is  not  obtained,  especially 
at  low  wind  speeds.  Even  with  a  good  first  guess  like  egn. 
(8.1.54),  we  found  that  simple  iteration  was  unreliable,  that  the 
Newton-secant  root  finding  procedure  was  needed  for  standard  cases, 
and  that  second  order  Aitken  acceleration  was  required  for  low  wind 
speeds  and  small  roughness  lengths. 

The  Newton-secant  root  finding  algorithm  utilizes  the  form, 

-<5,f  (x.) 

Xi+1     Xi  =  5i+l 


f(Xj)  -  f(XM)  (8.1.56) 

where  i  is  iteration  number.   Here  x  is  1/L,  and  f(x)  is  the 
difference  between  old  and  new  values  of  1/L.   The  object  of  the 
iteration  process  is  to  adjust  f(x)  toward  zero.   Unlike  the 
standard  Newton-Raphson  technique,  the  secant  method  does  not 
require  an  analytic  expression  for  the  first  derivative  of  f(x), 
which  in  our  case  is  not  obtainable.   The  Aitken  technique  is  a 
second  order  acceleration  method, 


Xi  +  3    ~    Xi    ~     5i  +  3    ~ 


XiXi  +  2  X  i+1 


xi  +  2    "    2xi  +  l    +    xi  (8.1.57) 


which  relies  on  the  Newton-secant  method  for  the  first  two 
iterations,  then  computes  the  rate  of  change  of  the  convergence 
from  the  first  two  iterations  and  uses  it  to  obtain  a  refined 
estimate  for  the  third  and  subsequent  iterations. 
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8.1.5   Objective  Analysis 

To  produce  each  of  the  n  flow  fields  for  the  n  towers,  the  mean 
wind  vector  from  each  tower  is  added  to  the  terrain  induced 
perturbation  field  from  LINCOM.  Then  at  each  grid  point,  the  Uj 
from  the  n  flow  fields  are  combined  according  to  the  algorithm, 

n 

S   U,/r^ 
i=l 

U(x,y)  =  ,  (8.1.58) 

n 

i=l 


where  r;  is  the  terrain  modified  distance  from  the  grid  point  to 
tower  i.   In  general  we  define  r  as, 

r  =  (a  +  (b  -  a)|sin  0|)R    ,  (8.1.59) 

Here  <p   is  the  tower-to-grid-point  angle  measured  from  north  (as  in 
meteorological  wind  angle) .  R  is  the  actual  tower-to-grid-point 
distance, 

a  =  V(A/B)    ,   and     b  =  1/a   . 

where  the  dominant  aspect  ratio  of  the  surrounding  terrain  height 
contours  is  given  by  A/B.  Thus,  rather  than  a  1/R  diminished 
circle  of  influence,  each  tower  has  a  1/r2  diminished  oval  of 
influence  which  corresponds  to  the  dominant  shape  and  orientation 
of  the  elevation  contours  surrounding  the  tower.  Of  course  for 
every  grid  point  the  r;  in  a  given  domain  need  be  computed  only 
once,  then  stored  for  future  use. 

This  method  of  combining  flow  fields  from  each  tower  ensures  that 
the  final  field  will  exactly  match  the  tower  data,  while  winds  at 
locations  between  towers  are  prescribed  by  a  combination  of  terrain 
modified  objective  analysis  and  neutral  boundary  layer  dynamics. 

In  LINCOM  1 . 1  we  noted  from  eqns.  (8.1.16  -  19)  that  a{   and  or2, 
the  vertical  decay  coefficients  for  the  terrain  induced  wind 
perturbation  fields,  are  independent  of  mean  wind  speed.   For 
neutral  flow  this  means  that  the  fundamental  characteristics  of  the 
perturbation  field  are  determined  by  the  terrain  and  are  thus 
domain  specific.   So  the  perturbation  field  may  be  expressed  as  a 
linear  combination  of  the  mean  wind  vector  and  the  two  orthogonal 
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component  solutions  for  the  perturbations  along  the  x  and  y  axes. 
Since  the  component  solutions  for  any  given  domain  may  be 
pre-stored,  this  allows  an  extremely  fast  "look-up  table"  mode  of 
running  LINCOM. 

Actually,  several  sets  of  orthogonal  solutions  may  be  stored  to 
account  for  different  inversion  heights.   Interpolation  and  a 
fitting  algorithm  are  then  used  to  determine  the  mean  wind  vector, 
which  when  combined  with  the  perturbation  fields,  will  best  fit  the 
available  tower  data.  The  best  fit  is  determined  by  minimizing  the 
differences  between  the  tower  and  LINCOM  winds.   Unlike  LINCOM  1.0 
this  new  procedure  in  LINCOM  1.1  relaxes  absolute  tower  anchoring 
in  order  to  ensure  mass  conservation  of  the  flow  field.    A 
reasonably  close  fit  to  the  tower  data  is  retained,  but  with  little 
potential  for  aliasing  as  described  in  sections  3.1  and  5.2.  This 
procedure  also  results  in  a  significant  gain  in  computational  speed 
when  many  towers  are  involved. 
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8.2    RIMPUFF  Model  Theory 


8.2.1   Introduction 

Since  distinctions  between  turbulence  and  mean  wind  are  not  well 
resolved  in  complex  terrain,  some  arbitrary  scale  such  as  the 
model's  grid  resolution  becomes  a  convenient  divider.   That  is,  we 
can  determine  the  mean  flow  down  to  the  grid  scale  as  a  direct 
deterministic  response  to  the  physical  forces.   However,  the  grid 
resolution  cannot  be  arbitrarily  fine  because  we  lack  fine  scale 
information  about  bottom  boundary  conditions,  especially  in  complex 
terrain,  and  also  because  computer  power  is  finite.   Moreover, 
sub-grid  scale  effects  must  be  parameterized  rather  than  simply 
neglected  because  the  non-linear  eddy  stress  dynamically  links  all 
scales  of  turbulence  in  an  energy  cascade.   Thus,  without  modeling 
the  dissipation  to  heat,  the  turbulence  kinetic  energy  would  simply 
grow  without  bound.   Luckily  for  our  investigation,  the  sub-grid 
fluid  motions  can  be  treated  statistically,  since  smaller  scales 
tend  to  behave  more  randomly,  even  in  complex  terrain. 

Hence,  RIMPUFF  lets  LINCOM  handle  the  bulk  transport,  while  it 
treats  the  instantaneous  diffusion  in  a  relative  frame  which 
follows  the  mean  flow.  Continuous  releases,  i.e.,  plumes,  are 
modeled  by  a  series  of  puffs  released  into  the  meandering  flow 
field.  RIMPUFF  computes  species  concentrations  by  summing 
contributions  from  all  puffs  for  each  grid  point  and  time  step. 
Depending  on  computer  capacity,  the  model  can  handle  several 
hundred  puffs  released  from  multiple  point  sources  anywhere  within 
the  domain.  Such  sources  can  have  individual  release  rates,  times 
and  heat  production  rates.  Plume  heights  depend  on  standard 
formulas  using  downwind  distance,  stability,  and  stability 
dependent,  vertical  profiles  of  windspeed.  The  inversion  is 
assumed  to  parallel  the  terrain  and  is  treated  as  a  material 
surface.  Both  the  inversion  and  source  heights  are  user  specified. 

Total  surface  reflection  of  puffs  is  normally  assumed  for  inert 
gases,  but  the  reflection  coefficient  is  adjustable  from  zero  to 
unity  for  materials  with  real  deposition  rates.  Dry  deposition  is 
calculated  from  source  depletion.  Dispersion  parameters  for  each 
puff  depend  on  items  such  as  pollutant  type,  atmospheric  stability, 
and  wind  speed.  Wet  deposition  depends  on  species  and  the  space/ 
time  distribution  of  the  rain  field. 


8.2.2   Governing  equations 

The  horizontal  sub-grid  diffusion  is  treated  by  linking  puff  growth 
to  the  time  averaged,  local  component  spectra  of  the  lateral 
velocities  (see  below) .   A  kinematic-statistical  approach  is 
employed  which  assumes  that  the  relative  displacement  of  fluid 
elements  proceeds  in  a  gaussian  manner. 
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That  is,  with  Q  referring  to  total  puff  mass  and  C(x,t)  to  the 
observed  concentration  field  in  space  and  time,  we  begin  by  posing 
the  mass  conservation  and  center  of  mass  location  as 

Q  -  jC(x,t) dx       ,    and  (8.2.1) 

c(t)  =  1/Q  Jx  C(x,t)3x    .  (8.2.2) 

The  center  of  mass  is  subject  to  transport  by  the  mean  flow 
determined  by  LINCOM  and  also  by  grid  scale  and  sub-grid  scale 
turbulence  not  included  in  LINCOM.  This  turbulence  "meander"  is 
estimated  by  a  stochastic  puff  advection  scheme  described  later. 
At  any  rate  the  center  of  mass  velocity  in  a  fixed  frame  can  be 
described  as  a  summation  of  the  product  of  the  puff  mass  at  each 
location  with  its  respective  velocity, 

Vcm(t)   =   1/Q  J  u(x,t)  C(x,t)  ax     ,  (8.2.3) 


where  u(x,t)  represents  the  fluid  velocity  field.   We  then  denote 
both  the  average  of  all  particles  within  a  puff  and  the  ensemble 
average  by  "<  >"  and  fluctuations  by  "  •  ".   Thus,  if  C(x,t)  = 
<  C(x,t)  >  +  C'(x,t),  then  1/Q<  C(x,t)  >dx  denotes  the  probability 
of  finding  a  marked  particle  at  the  point  (x,t)  and 

<  x2  >  =  1/Q  J  x2<  C(x,t)  >dx  (8.2.4) 

denotes  the  ensemble  averaged,  fixed  frame,  absolute  diffusion. 
The  latter  may  be  expressed  as  the  sum  of  meander  and  relative 
diffusion,  or 

<  x2  (t)  >  =  <  y2  (t)  >  +  <  c2  (t)  >    ,  (8.2.5) 

where  y  =  x  -  c  defines  a  coordinate  frame  moving  with  the  center 
of  mass  velocity. 

Now,  using  t  to  time  integrate  along  a  particle's  path  with 
respect  to  the  center  of  mass,  the  puff  growth  rate  becomes, 

t 
do2/dt   =      2\    <    v(t)  v(t-r)  >  dT 
o 
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=   <  v2(t)  >\    r(t,r)dr   =   <  v2  (t)  >  tr(t) 
o 


(8.2.6) 


Here,  a  is  puff  size  (standard  deviation  of  particle  displacement 
relative  to  center  of  mass)  and  <  v2  (t)  >  is  the  corresponding 
relative  velocity  variance.  The  relative  correlation  function, 
r(t,T),  depends  on  both  t  and  t.  This  indicates  a  non-stationary 
process,  unlike  Taylor's  single  particle  diffusion.  Letting  v  =  u 
-  Vcm  be  the  moving  frame  particle  velocity,  we  can  show  that 


do*/dt  =   2  J  [Rabs(0  -  Rcm(t,T)]  3t 
o 


(8.2.7) 


Here, 


Rabs(T)    =  <  u(t)u(t-T)  >    ,     and 


(8.2.8a) 


Rcm(t,T)  =  <  Vcm(t)Vcm(t-r))  > 


(8.2.8b) 


are  the  appropriate  auto-covariance  functions  for  absolute  single 
particle  diffusion,  and  from  the  center-of-mass-velocity  viewpoint, 


respectively 
for  R  . 


Since  R  b<i  is  known,  we  can  focus  on  obtaining  a  model 


Using  (8.2.3),  R  (t,r)  can  be  rewritten  as 


H  <  u(x' ,t)u(x",t-T)  C(x' ,t)C(x",t-T)  >dx'dx" 


(8.2.9) 


From  this  point  the  puff  can  be  modeled  as  a  set  of  passive  marked 
particles  with  instantaneous  concentration,  C(x,t),  moving  within 
a  large  but  numerable  set  of  unmarked  fluid  elements,  with  Eulerian 
velocity  field  u(x,t).  Any  relative  displacements,  y(t),  can  also 
be  assumed  to  obey  the  same  un-correlated  gaussian  statistics  (see 
fig.  8.2.1).  That  is,  any  two  particles  viewed  in  the  moving  frame 
will  appear  to  disperse  in  an  uncorrelated  fashion,  such  that 


<Ay   Ay  >  = 
i    J 


a2  (t)  -  a2  (t-r)    for  i  =  j 
0  for  i  ^  j 


(8.2.10) 
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Physically,  this  means  that  the  characteristic  scale  of  motion  in 
the  moving  frame  must  be  much  shorter  than  the  puff  size.  However, 
these  scales  may  overlap  in  the  atmosphere,  implying  that  only  the 
ensemble  average  of  individual  puff  releases  will  be  approximately 
gaussian. 

Hence,  a  gaussian  puff  of  initial  size,  a(0),  will  remain  a 
gaussian  puff  of  size,  a(t4),  at  a  later  time,  (t;)  ,  so  that 

Gff(t)(y,t)  =  J  Gs(y-y0,t)  Go(t_r)(y0,t-T)3y0    .  (8.2.11) 


We  also  see  that  the  spread,  <  Ay2j(t,T)  >,   for  the  transition 
probability,  Gs,  given  in  eqn.  (8.2.10)  satisfies  eqn.  (8.2.11). 

Substituting  the  instantaneous  concentrations  from  eqn.  (8.2.9) 
with  the  corresponding  gaussians  now  leads  to 


Rcmft'T)  =  H<  u(y\t)u(y",t-T)  >G0(t)(y\t)Go(tr)(y",t-T)ay'dy", 

(8.2.12) 

where  the  moving  frame  fluid  velocity,  u,  relates  to  the  fixed 
frame  velocity  by  u(y,t)  =  u(y+c,t). 

We  must  also  obtain  the  two-point,  two-time  (fixed  frame)  velocity 
covariance,  <  u (y ' , t) u(y" , t-T )  >.  For  two  particles,  i  and  j,  in 
a  homogeneously  turbulent  field,  the  velocity  covariance  function, 
in  a  fixed  frame,  Rahs  (£-,  t )  ,  will  depend  only  upon  the  temporal  and 
spatial  separations,  i  and  £»,  between  the  particles.   That  is, 


RabS(Sij'T)  =  <  u(x,(t))  u(x,(t-r)  -  £ij)  >    .  (8.2.13) 

Thus,  the  moving  frame  velocity  covariance  in  eqn.  (8.2.12)  is 
written  in  terms  of  the  absolute  velocity  covariance  function  and 
the  transition  probability  from  eqn.  (8.2.11)  as 

<  u(y\t)u(y»,t-r)  >  =  j  Gs(y»-y"-£,t)RabsU,T)d£    .   (8.2.14) 

Returning  this  to  eqn.  (8.2.12)  and  integrating  over  y'  and  y"  will 
reformulate  the  center-of-mass  covariance  as 
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Rcm(t,T)  -  Rabs(^T)/[2a(t)V/7r  e-^/4ff,W]3|   .  (8.2.15) 

Thus,  in  terms  of  the  fixed  frame  covariances,  the  horizontal 
rate  equation  for  puff  growth  becomes, 

t 

3o*/dt    =  2  J    (Rabs(0,T)  - 

o 


Rabs(£,T)/[  2Vtt  tJ(t)  e«,/4^]  d^T    . 


(8.2.16) 


For  ease  of  use,  we  employ  the  Fourier  transform, 

Rahs(£,T)  =  S(K,o>)    e^+^dco    ,  (8.2.17) 

to  render  eqn. (8.2.16)  as, 


do*/8t   =    2tJJ    S(k,w)    sin(wt)/wt    (1    -    e"K?a?(,))     dcodx 


(8.2. 18) 


In  this  case  the  spatial  filter,  1  -  e"*o(t),  has  been  inserted  to 
remove  wave  numbers  smaller  than  1/a,  while  the  low  pass  temporal 
frequency  filter,  sin  (cot)  /cot,  does  likewise  beyond  co  =  1/t.  The 
shaded  area  in  Fig.  8.2.2  shows  the  part  of  S(k,co)  that  effectively 
contributes  to  the  growth  rate. 

For  puffs  much  larger  than  the  turbulence  length  scale,  eqn. 
(8.2.18)  reduces  to  Taylor's  formula  for  single  particle  diffusion. 
Hence,  behavioral  differences  between  a  puff  and  a  single  particle 
are  closely  related  to  the  spatial  correlation  of  the  turbulence. 
That  is,  for  classically  homogeneous  turbulence  with  a  Kolmogorov 
inertial  subrange,  the  spectral  index,  p  =  -5/3.  Puff  growth  will 
then  follow  a  standard  t3/2  prediction.  However,  in  a  stably 
stratified  atmosphere,  p  may  be  =  -3,  implying  exponential  growth. 
In  either  case,  we  assume  that  the  sub-grid  diffusion  is  locally 
homogeneous.  If  so,  we  can  scale  the  puff's  instantaneous  growth 
rate  over  intermediate  ranges  (<  500m)  with  the  local  lateral 
turbulence  intensity,  ot,  such  that  do/dt  ~  0.3ajU.  We  obtained 
turbulence  intensities  for  Vandenberg  by  modeling  Eulerian  wave 
lengths  with  the  similarity  formula, 
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a-      =  oi0(  1  +  B/ut)1/3        ,  (8.2.19) 


where  o-0  refers  to  limiting  values  at  large  times  (Skupniewicz ,  et 
al,  1989).  That  is,  even  within  complex  terrain,  the  turbulent 
energy  characterized  by  Vandenberg *s  tower  wind  spectra  tends  to 
plateau  at  larger  scales.  Eguation  (8.2.19)  tends  to  apply  over  a 
wider  range  of  scales  than  simple  power  law  estimates.  Power  law 
fits  to  Vandenberg  data  are  also  found  to  be  rather  windspeed 
dependent.  Note  that  the  constant,  B,  depends  upon  measurement 
height  and  stability  and  is  proportional  to  the  dominant  eddy  size 
at  a  given  height.  Thus,  values  for  ai0  and  B  were  determined  by 
multi-variate  regression  for  each  tower  for  each  of  three  periods: 
pre-dawn,  noon,  and  dusk,  indicative  of  stable,  unstable,  and 
neutral  conditions.  Values  were  also  classed  by  wind  direction. 
The  longitudinal  and  lateral  values  determined  for  ai0  ranged 
typically  from  700  -  2,000m  and  200  -  1,000m,  respectively.  As  in 
flat  terrain,  the  longitudinal  component  clearly  dominates. 
However,  our  values  were  consistently  larger,  suggesting  that  flow 
distortions  due  to  complex  terrain  feed  energy  directly  into  larger 
horizontal  scales  of  motion. 


8.2.3   Puff  splitting 

Real  diffusion  in  complex  terrain  also  often  displays  plume 
bifurcation  and/or  vertical  shear  due  to  channeling,  slope  flows, 
or  inversions.  To  simulate  such  decoupling  we  allow  each  puff  to 
split  as  often  as  necessary.  The  daughter  puffs  are  arranged  such 
that  the  original  center  concentration,  radially  integrated  second 
moments  of  the  concentration  distribution,  and  mass  are  all 
conserved . 

These  constraints  may  be  summarized  as, 


2 

iS  =   a2pl     ,                           (8.2.20a) 
5 

ii     C5(0,0)  =   0,(0,0)     ,                         (8.2.20b) 

iii    o?5  =   (l/2)crpl      ,    and                 (8.2.20c) 

iv  mass  conservation. 
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For  the  Vandenberg  simulations  the  initial  puff  diameter  was  set  to 
100m  (  one  standard  deviation  in  concentration) .   When  the  puff 
attains  a  diameter  equal  to  the  500m  grid  spacing,  it  splits 
(pentafurcates)  horizontally  into  five  new  gaussian  puffs.   The 
daughter   puff   diameters   are   initially   set   to   250m.     The 
conservation  relations  then  require  that  the  four  satellite  puffs 
be  centered  at  0.89  ctv  from  the  origin  and  that  the  satellite  and 
center  puffs  each  carry  23.53  and  5.88%  of  the  total  matter, 
respectively  (fig. 5). 

In  the  case  of  layers  decoupled  vertically  by  shear,  the  original 
puff  is  allowed  to  trifurcate  into  three  daughter  puffs  which  are 
centered  along  the  original  puff's  vertical  axis.  To  simulate  a 
gaussian  radial  concentration  profile,  the  conservation  rules 
require  that  the  two  satellite  puffs  be  centered  at  +/-1.19av,  with 
each  new  puff  carrying  26.58  %  of  the  initial  pollutant  mass. 
The  new  center  puff  retains  the  remaining  48.84  %.  Again,  the  new 
puffs  are  born  with  half  the  diameter  of  the  original  puff. 

The  purpose  of  this  splitting  scheme  is  to  allow  each  new  puff  its 
own  individual  growth  rate  and  direction,  as  set  by  the  mean  flow 
and  turbulence  intensities  local  to  it.   We  assume  in  this 
procedure  that  the  sub-grid  scale  shearing  is  locally  homogeneous 
so  that  filtered  turbulence  estimates  can  be  used  in  the 
parameterization.   The  puff-splitting  scheme  can  be  repeated  to  a 
practical  limit  of  several  hundred  puff  progeny.  In  this  way  the 
plume  bifurcation  caused  by  terrain  features  can  be  modeled  in 
detail.   Such  bifurcation  was  seen  extensively  during  LVDE  (see 
LVDE  Data  Report,  1990) 

8.2.4    Stochastic  puf f-advection 

The  mean  flow  field  provided  by  LINCOM  should  be  updated  every  5  to 
10  minutes  to  include  temporal  variability  and  allow  for  puff 
meander  in  RIMPUFF  due  to  non-stationary  winds.  But  if  measured 
winds  are  available  only  as  one  hour  means,  only  horizontal  shear 
in  the  static  wind  field  is  left  to  induce  any  meander  in  the  puff 
center-of-mass  advection  velocity,  Vcm.  Such  was  the  case  for  the 
eight  Mt.  Iron  diffusion  experiments.  To  re-introduce  temporal 
variability  in  the  wind  field,  we  have  added  an  auto-regressive 
(Langevin  equation  based)  advection  component  to  the  two  horizontal 
(u,  v)  wind  components  of  the  mean  flow  model.  This  is  akin  to 
standard  Langevin  random  forcing  in  Monte  Carlo  dispersion  but 
applied  to  puffs  rather  than  particles. 

We  estimated  supra-grid  scale  (>  500  m)  wind  variance,  <VcmVcm>,  by 
low-pass  filtering  the  horizontal  wind  energy  spectrum  below  wave 
number  k  =  2a/500  m.  Since  we  already  use  the  sub-grid  scale  (< 
500  m)  variance  for  puff  growth,  we  now  account  for  turbulence  at 
all  scales.   TL   ff  =  LE/u  also  provides  a  Lagrangian  time  scale  for 
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horizontal  puff  motion,   where  LE,  an  Eulerian  length  scale,  is 
derived  from  velocity  spectra  measurements  (see  Handbook) . 

With  these  we  can  estimate  stochastic  contributions  to  Vcm   from 

Vcm(t  +  At)  =  PVcm(t)  +  V        .  (8.2.21) 

where  the  Lagrangian  correlation  coefficient  for  the  center-of-mass 
velocity  is  given  by 

p    =    exp(-At/TL)     ,  (8.2.22) 

and  r\    is  white  gaussian  noise  having  a  variance,  (1  -  P2)  <VcmVcm>. 

By  comparing  cases  with  and  without,  we  found  that  the  stochastic 
component  is  clearly  visible  during  individual  simulations. 
However,  its  contribution  to  the  overall  plume  widths  is  small. 
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TIME 


Fig.  8.2.1  The  moving  frame  trajectory,  y;,  of  a  marked  fluid 
particle,  i,  that  at  time  t-T,  holds  the  position,  yj(t-T).  The 
guantity,  Ay,  s  y^t)  -  Vj(t-T)  ,  and  its  gaussian  distribution 
function,  Gs,  are  shown  at  time  t. 
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Fig.  8.2.2  Schematic  isopleth  plot  of  spectrum,  S(k,oj).  Its 
maximum  value  is  at  (k,fo)  =  (0,0)  from  where  the  function  decreases 
monotonically  through  the  levels  I,  II,  and  III.  The  shaded  region 
is  the  effective  domain  of  the  model  defined  by  vertical  line,  b  = 
t"1,  (using  the  low-pass  filter,  sin (fot) /fot)  and  the  horizontal 
line,  k  =  o'x    (using  the  high  pass  filter,  1  -  exp(-k2a2)  ). 
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Pen  I  a  Puff 


Fig.  8.2.3  Puff  pentaf urcation:  the  original  single  puff  of  size 
a  divides  into  five  new  puffs,  each  of  size  \  o,  conserving  total 
mass,  and  peak  and  second  moments  of  the  concentration  distribution 
around  the  cloud's  center  of  mass  point  (moment  of  inertia). 
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8.3    Comparison  Procedures 

The  following  describes  how  the  LINCOM/RIMPUFF  model  runs  were 
compared  with  the  digitized  Mt.  Iron  data  for  cases  28,  31,  48,  55, 
87,  90,  91  and  110. 

First,  the  LINCOM/RIMPUFF  model  runs  and  Mt.  Iron  case  studies  were 
plotted  with  the  SURFER  graphics  package  so  the  plumes  could  be 
visually  inspected.  By  overlaying  a  model  run  with  a  case  study  it 
was  evident  for  some  cases  that  the  source  points  were  not 
collocated.  The  Mt.  Iron  grids  were  moved  by  an  appropriate  number 
of  grid  points  in  the  x  and  y  directions  to  align  the  source  point 
of  the  Mt.  Iron  case  study  and  the  LINCOM/RIMPUFF  model  run. 

Headers  at  the  top  of  the  SURFER  formatted  *.grd  files  define  the 
latitude/ longitude  coordinates  in  UTMS  for  the  four  corners  of  the 
LINCOM/RIMPUFF  and  Mt .  Iron  dose  files.  For  all  of  the  LINCOM/ 
RIMPUFF  runs,  the  headers  did  not  correctly  locate  the  source 
points  relative  to  a  map  of  Vandenberg.  By  plotting  some  of  the 
roads  at  Vandenberg  along  with  the  location  of  VIP-1,  the  proper 
coordinates  for  the  surfer  header  files  for  each  of  the  model  runs 
were  determined.  These  roads  were  plotted  by  picking  off  latitude 
and  longitude  coordinates  for  points  on  the  roads  using  a  300'  (.3 
arc  seconds)  resolution  map  of  Vandenberg.  Program  utm.for  was 
written  to  convert  the  longitude/latitude  values  to  UTMS 
coordinates . 

For  cases  28,  48,  87,  90  and  110  the  plotted  roads  (which  were  also 
sampling  lines  for  the  Mt.  Iron  data  set)  were  used  to  define  plume 
cutoff  lines.  Since  the  plumes  did  not  reach  Santa  Ynez  Rd .  in 
cases  28  and  110,  these  plumes  were  truncated  at  Arguello  Rd.  The 
plumes  were  truncated  at  Santa  Ynez  Rd.  for  cases  87  and  90  because 
there  was  enough  data  to  define  these  plumes  at  Santa  Ynez.  For 
case  48,  the  plume  was  cut  off  at  Honda  Ridge. 

The  following  procedure  was  used  for  this  evaluation.  1)  TEST,  an 
ASCII  file  was  created  to  hold  the  input  and  output  names  of  the 
LINCOM/RIMPUFF  and  Mt .  Iron  grid  files.  This  file  is  comprised  of 
four  file  names  followed  by  a  blank  line  for  each  comparison.  For 
example,  to  compare  a  digitized  data  grid  with  a  model  simulation 
grid  for  case  28,  TEST  would  look  like: 

mi28.grd  #  input  file  name  of  Mt.  Iron  data  grid 

case28.grd  #  input  file  name  of  LINCOM/RIMPUFF  simulation  grid 

p28.grd  #  output  file  name  of  Mt.  Iron  grid 

c28.grd  #  output  file  name  of  LINCOM/RIMPUFF  grid 

2)  The  program,  drawsurfl.f,  was  run  to  collocate  the  release 
points,  find  the  merit  scores  where  the  hand-drawn  isopleths  were 
left  untruncated  (cases  31,  55  and  91),  and  to  create  an  output 
file  to  contour  with  TOPO. 
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3)  The  header  files  in  all  of  the  LINCOM  model  files  were  changed 
so  that  the  x  and  y  are  in  UTMS  coordinates: 

4)  Cases  31,  55  and  91  were  then  complete  and  ready  to  be  plotted. 

5)  Program  mt3.for  was  run  for  the  truncation  cases:  28,  87,  90  and 
110.  All  concentrations  beyond  the  last  sampling  line  where 
appreciable  doses  were  measured  were  set  to  zero. 

6)  mt4 . f or  was  run  for  case  48  which  was  truncated  at  Honda  Ridge. 

7)  Program  drawsrf2.f  was  run  on  cases  28,  48,  87,  90  and  110. 
This  program  is  actually  the  same  as  drawsrfl.f  except  that  the 
format  statements  to  read  the  input  grids  are  different  and  also 
the  source  point  adjustment  is  not  re-evaluated. 

8)  Cases  28,  48,  87,  90  and  110  were  then  complete  and  plotted. 
The  "roads. bin"  file  was  used  to  overlay  some  roads  on  the  plots. 
Some  of  these  roads  correspond  to  sampling  lines. 
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